Digit DP 5 trạng thái. Ta gọi hàm digit DP của mình là ~\texttt{dp(pos, geqL, leqR, bound, diff_mod)}~ với

• ~\texttt{pos}~: Vị trí chữ số đang xét trong hệ nhị phân của 2 số ~a,b~
• ~\texttt{geqL}~: ~a>l~ hay chưa
• ~\texttt{leqR}~: ~b < r~ hay chưa
• ~\texttt{bound}~: ~a < b~ hay chưa
• ~\texttt{diff_mod}~: Chênh lệch số dư khi chia cho ~k~ của ~a~ và ~b~ hiện tại.

Việc còn lại là tiến hành digit DP chuyển trạng thái và xét các trường hợp ~0/1~ tương ứng của các cặp số ~(a,b)~.

Code mẫu

#include <bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>

#define fore(i, a, b) for (int i = (a), _b = (b); i < _b; ++i)
#define fort(i, a, b) for (int i = (a), _b = (b); i <= _b; ++i)
#define fordt(i, a, b) for (int i = (a), _b = (b); i >= _b; --i)
#define forde(i, a, b) for (int i = (a), _b = (b); i > _b; --i)
#define trav(a, x) for (auto& a : x)
#define sz(x) ((int)(x).size())
#define all(x) (x).begin(), (x).end()

using namespace std;
using namespace __gnu_pbds;

template<class A, class B> bool maxi(A &a, const B &b) {return (a < b) ? (a = b, true):false;};
template<class A, class B> bool mini(A &a, const B &b) {return (a > b) ? (a = b, true):false;};

typedef int64_t i64;
typedef pair<int,int> _ii;
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> ordered_set;
constexpr int maxn=1e4+7;
constexpr i64 oo=1e9+7;
i64 dp[65][2][2][2][maxn], l, r, pw[65];
int k;

i64 calc(int pos, int geqL, int leqR, int bound, i64 diff_mod) {
    if(dp[pos][geqL][leqR][bound][diff_mod]!=-1)
        return dp[pos][geqL][leqR][bound][diff_mod];

    if(!pos)
        return dp[pos][geqL][leqR][bound][diff_mod] = (diff_mod == 0) & bound;

    i64 ans=0;
    int limL=(l >> (pos - 1))&1;
    int limR=(r >> (pos - 1))&1;
    fort(dA, 0, 1) fort(dB, 0, 1) {
        if(!bound && dA > dB) continue;
        if(!geqL && dA < limL) continue;
        if(!leqR && dB > limR) continue;

        int geqL_new = geqL;
        if(dA > limL) geqL_new = 1;
        int leqR_new = leqR;
        if(dB < limR) leqR_new = 1;
        int bound_new = bound;
        if(dA < dB) bound_new = 1;

        int bit_and = (dA & dB) & 1, bit_xor = (dA ^ dB) & 1;
        i64 diff_mod_new = diff_mod;
        if(bit_and) diff_mod_new = (diff_mod_new + pw[pos - 1]) % k;
        if(bit_xor) diff_mod_new = (diff_mod_new - pw[pos - 1] + k * oo) % k;

        ans = (ans + calc(pos - 1, geqL_new, leqR_new, bound_new, diff_mod_new)) % oo;
    }

    // cout << ans << endl;
    return dp[pos][geqL][leqR][bound][diff_mod] = ans;
}

i64 brute() {
    i64 ans=0;
    fort(i, l, r) fort(j, i + 1, r) {
        i64 a=i^j, b=i&j;
        if((a % k) == (b % k)) {
            cout << i << ' ' << j << ": " << a << ' ' << b << '\n';
            ++ans;
        }
    }

    // cout << ans << '\n';
    return ans;
}

void solve() {
    cin >> l >> r >> k;
    assert((l <= r) && (1 <= k) && (k <= 1e4));

    pw[0] = 1;
    fort(i, 1, 64) pw[i] = (pw[i - 1] << 1) % k;

    memset(dp, -1, sizeof(dp));
    i64 ans = calc(62, false, false, false, 0);
    // assert(ans==brute());
    // cout << brute() << ' ' << ans << '\n';
    cout << ans << '\n';
}

#define NAME "test."
int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(0); cout.tie(0);

    if(fopen(NAME"inp", "r")) {
        freopen(NAME"inp", "r", stdin);
        freopen(NAME"out", "w", stdout);
    }

    int t=1;
    while(t--) solve();
}