SP26017 GCDMAT

求：

$$\sum_{i=i_1}^{i_2}\sum_{j=j_1}^{j_2}\gcd(i,j) \pmod{10^9+7}$$

$T$ 组数据。

其中，$T \leq 50000,1 \leq i2,j2\leq 10^6$。

首先设 $S(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\gcd(i,j)$。

则 $ans=S(i_2,j_2)+S(i_1,j_1)-S(i_2,j_1)-S(i_1,j_2)$。

考虑化简 $S(n,m)$，有：

$$\large \begin{aligned} S(n,m)&=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\gcd(i,j)\\ &=\sum_{d=1}^{n}d\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=d]\\ &=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)=1] \\ &=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{p|\gcd(i,j)}^{}\mu(p)\\ &=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)[p|i][p|j] \\ &=\sum_{d=1}^{n}d\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[p|i][p|j] \\ &=\sum_{d=1}^{n}d\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\lfloor\frac{n}{pd}\rfloor\lfloor\frac{m}{pd}\rfloor \\ &=\sum_{d=1}^{n}d\sum_{p=1}^{\lfloor\frac{n}{d}\rfloor}\mu(p)\lfloor\frac{n}{pd}\rfloor\lfloor\frac{m}{pd}\rfloor \\ &=\sum_{T=1}^{n}\sum_{p|T}\mu(p)\frac{T}{p}\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor & 设\ T=dp\\ &=\sum_{T=1}^{n}\varphi(T)\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor \\ \end{aligned}$$

这是 $\mathcal{O}(T\sqrt{n})$ 的算法，且常数巨大。

考虑优化，发现：

$$\large \begin{aligned} ans=\sum_{i=1}^{n}\varphi(i)(\lfloor\frac{i_2}{i}\rfloor-\lfloor\frac{i_1}{i}\rfloor)(\lfloor\frac{j_2}{i}\rfloor-\lfloor\frac{j_1}{i}\rfloor)&&n=\min(i_2,j_2) \end{aligned}$$

但是，还是过不了。

这里有一个好方法：预处理出 $1\sim N$ 的倒数，然后把除法改成乘法，优化常数，可过。

但还可以优化。

考虑根号分治。

我们发现这实际上是一个预处理和查询的平衡，因为我们发现，$l$ 和 $r$ 越大，相同一段长度就越大，暴力的复杂度就相对越劣（意思就是用整除分块处理越快），那么我们可以考虑在 $l$ 和 $r$ 较小的时候暴力，否则预处理。

详见代码。

#pragma GCC optimize("Ofast")

#include <cstdio>

#include <cmath>

#define int long long

const int p = 1e9 + 7;

namespace Fread
{
    const int SIZE = 1 << 21;
    char buf[SIZE], *S, *T;
    inline char getchar()
    {
        if (S == T)
        {
            T = (S = buf) + fread(buf, 1, SIZE, stdin);
            if (S == T)
                return '\n';
        }
        return *S++;
    }
    inline int read()
    {
        int x = 0, f = 1;
        char c = getchar();
        while (c < '0' || c > '9')
        {
            if (c == '-')
                f = -1;
            c = getchar();
        }
        while (c >= '0' && c <= '9')
        {
            x = (x << 1) + (x << 3) + (c ^ 48);
            c = getchar();
        }
        return x * f;
    }
}

namespace Fwrite
{
    const int SIZE = 1 << 21;
    char buf[SIZE], *S = buf, *T = buf + SIZE;
    inline void flush()
    {
        fwrite(buf, 1, S - buf, stdout);
        S = buf;
    }
    inline void putchar(char c)
    {
        *S++ = c;
        if (S == T)
            flush();
    }
    struct NTR
    {
        ~NTR() { flush(); }
    } ztr;
    inline void write(int n, char c)
    {
        int tot = 0, a[101];
        while (n)
            a[++tot] = n % 10, n /= 10;
        if (!tot)
            putchar('0');
        else
            while (tot)
                putchar(a[tot--] + 48);
        putchar(c);
    }
}

int min(int n, int m)
{
    return n < m ? n : m;
}

int max(int n, int m)
{
    return n > m ? n : m;
}

int pri[10000001], sphi[10000001], phi[10000001], tot, N, mid;

double inv[10000001];

bool vis[1000001];

int opt;

int ge(int x, int y)
{
    if ((int)(x * inv[y]) != 0)
        return (int)(x * inv[(int)(x * inv[y])]);
    return opt;
}

void solve()
{
    sphi[1] = phi[1] = 1;
    for (int i = 2; i <= N; ++i)
    {
        if (!vis[i])
            pri[++tot] = i, sphi[i] = i - 1;
        for (int j = 1; j <= tot && pri[j] * i <= N; ++j)
        {
            vis[pri[j] * i] = 1;
            if (!(i % pri[j]))
            {
                sphi[pri[j] * i] = 1ll * sphi[i] * pri[j];
                break;
            }
            sphi[pri[j] * i] = 1ll * sphi[i] * (pri[j] - 1);
        }
        phi[i] = sphi[i];
    }
    for (int i = 1; i <= N; ++i)
        sphi[i] += sphi[i - 1], inv[i] = (1.00000000 / i) * (1 + 1e-15);
}

signed main()
{
    int t = Fread::read();
    N = max(Fread::read(), Fread::read());
    solve();
    while (t--)
    {
        int x1 = Fread::read() - 1, y1 = Fread::read() - 1, x2 = Fread::read(), y2 = Fread::read(), ans = 0;
        if (x2 > y2)
        {
            x1 ^= y1 ^= x1 ^= y1;
            x2 ^= y2 ^= x2 ^= y2;
        }
        opt = x2;
        mid = sqrt(y2 * 7);
        for (int i = 1; i <= mid; i++)
            ans += phi[i] * ((int)(x2 * inv[i]) - (int)(x1 * inv[i])) * ((int)(y2 * inv[i]) - (int)(y1 * inv[i]));
        ans %= p;
        for (int l = mid + 1, r; l <= x2; l = r + 1)
        {
            r = min(min(ge(x1, l), ge(y1, l)), min(ge(x2, l), ge(y2, l)));
            ans += 1ll * (sphi[r] - sphi[l - 1]) * ((int)(x2 * inv[l]) - (int)(x1 * inv[l])) * ((int)(y2 * inv[l]) - (int)(y1 * inv[l]));
            ans %= p;
        }
        Fwrite::write(ans, '\n');
    }
    return 0;
}