Gọi ~f(x, y, k)~ là tổng các tích ~i \times j~ của các cặp ~i, j~ thỏa ~i \leq x, j \leq y~ và ~gcd(i, j) = k~.

Vậy đáp án cho bài toán là ~f(b, d, k) - f(a - 1, d, k) - f(b, c - 1, k) + f(a - 1, c - 1, k)~

Ta có ~f(x, y, k) = f(x / k, y / k, 1) \times k^2~, bài toán chuyển về tính ~f(x, y, 1)~.

Công thức tính ~f(x, y, 1)~ như sau:

~f(x, y, 1) = \displaystyle\sum_{k = 1}^{min(x, y)}{ g(x / k, y / k) \times mo[k] \times k^2 } ~, trong đó, ~mo~ là hàm mobius, ~g(a, b)~ là tổng các tích ~i \times j~ của các cặp ~i, j~ mà ~i \leq a, j \leq b~, ta có ~g(a, b) = (a \times (a + 1) / 2) \times (b \times (b + 1) / 2)~.

Nhận xét rằng, có thể chia đoạn ~[1;min(x, y)]~ thành khoảng ~\sqrt{min(x, y)}~ đoạn ~[l; r]~ mà ~x / k~ không đổi và ~y / k~ không đổi với mọi ~k~ thuộc đoạn ~[l;r]~.

Với mỗi đoạn ~[l;r]~ như vậy, ta có ~\displaystyle\sum_{k = l}^{r}{ g(x / k, y / k) \times mo[k] \times k^2 = g(x / l, y / l) \times \displaystyle\sum_{k = l}^{r}{mo[k] \times k^2} } ~, phần ~(\displaystyle\sum_{k = l}^{r}{mo[k] \times k^2})~ có thể được tính trong ~O(1)~ bằng cách tạo trước mảng tổng tiền tố của hàm ~mo[k] \times k^2 ~.

Vậy mỗi đoạn ~[l; r]~ chỉ cần ~O(1)~ để tính toán, độ phức tạp cho mỗi bộ dữ liệu là ~O(\sqrt{min(b/k, d/k)})~ hay ~O(\sqrt{min(b, d)})~, tổng độ phức tạp là ~O(q\sqrt{b})~.

Code mẫu

#include <bits/stdc++.h>

#define fi first
#define se second
#define For(i, a, b) for (int i=a;i<=b;++i)
#define Ford(i, a, b) for(int i=a;i>=b;--i)

using namespace std;

typedef long long ll;
typedef pair<int, int> ii;
typedef vector<int> vi;
typedef vector<ll> vll;
typedef vector<ii> vii;

const int N = 2 * 1e5 + 5;
const ll mod = 1e9 + 7;

ll p[N], mu[N], m2[N];

int mobius(int x) {
    int c = 0;
    while (x > 1) {
        if (x / p[x] % p[x] == 0) return 0;
        c++;
        x /= p[x];
    }
    return (c & 1 ? -1 : 1);
}
void prepare(int n) {
    For(i, 1, n) p[i] = i;
    for (int i = 2; i * i <= n; ++i) {
        if (p[i] == i) {
            for (ll j = i * i; j <= n; j += i) {
                p[j] = i;
            }
        }
    }

    For(i, 1, n) m2[i] = (m2[i - 1] + 1LL * i * i %mod * (mobius(i) + mod)) % mod;
}

ll calc(ll x, ll y) {
    return x * (x + 1) / 2 % mod * (y * (y + 1) / 2 % mod) % mod;
}
ll f(int x, int y) {
    if (x > y) swap(x, y);
    ll ans = 0;

    int l = 1;
    while (l <= x) {
        int r = min(x / (x / l) + 1, y / (y / l) + 1);

        (ans += calc(x / l, y / l) * (m2[r - 1] - m2[l - 1] + mod)) %= mod;
        l = r;
    }
    return ans;
}

void solve() {
    int a, b, c, d, g;
    cin >> a >> b >> c >> d >> g;
    if (g < 1) cout << "0\n";
    else
        cout << (f(b / g, d / g) - f((a - 1) / g, d / g) - f(b / g, (c - 1) / g) + f((a - 1) / g, (c - 1) / g) + mod + mod) * g % mod * g % mod << "\n";
}

int main() {
//  freopen(".inp","r",stdin);
//  freopen(".out","w",stdout);

    ios_base::sync_with_stdio(false);
    cin.tie(NULL);cout.tie(NULL);

    prepare(100000);

    int t; cin >> t;
    while (t--) solve();

    return 0;
}