BZOJ 4018 小Q的幻想之乡

题目：http://www.lydsy.com/JudgeOnline/problem.php?id=4018

题意：
$T$组询问，每组询问给定$N$和$M$，求$\sum_{i=1}^{N}{\sum_{j = 1}^{M}{\frac{|i-j|}{\gcd(i,j)}}}$分别模$10^9+7$和$10^9+9$的值。
$T\le1000, N, M \le2*10^6$。

题解：
考虑将$\sum_{i=1}^{N}{\sum_{j = 1}^{M}{\frac{|i-j|}{\gcd(i,j)}}}$化简。
首先引入一个记号，令$[x]$表示：若$x$为真，则$[x]$为$1$；否则$[x]$为$0$，其中$x$是一个布尔表达式。

$$\sum_{i=1}^{N}{\sum_{j = 1}^{M}{\frac{|i-j|}{\gcd(i,j)}}}\= \sum_{i=1}^{N}{\sum_{j = 1}^{M}{\sum_{d}{\frac{|i-j|}{\gcd(i,j)}[d=\gcd(i,j)]}}}\= \sum_{d}{\sum_{i=1}^{\left \lfloor \frac{N}{d} \right \rfloor}{\sum_{j = 1}^{\left \lfloor \frac{M}{d} \right \rfloor}{|i-j|[\gcd(i,j)=1]}}}\= \sum_{d}{\sum_{i=1}^{\left \lfloor \frac{N}{d} \right \rfloor}{\sum_{j = 1}^{\left \lfloor \frac{M}{d} \right \rfloor}{|i-j|\sum_{d‘|\gcd(i,j)}{\mu(d‘)}}}}\= \sum_{d}{\sum_{d‘}{d‘\mu(d‘)\sum_{i=1}^{\left \lfloor \frac{N}{dd‘} \right \rfloor}{\sum_{j = 1}^{\left \lfloor \frac{M}{dd‘} \right \rfloor}{|i-j|}}}}\= \sum_{d}{\sum_{d‘|d}{d‘\mu(d‘)\sum_{i=1}^{\left \lfloor \frac{N}{d} \right \rfloor}{\sum_{j = 1}^{\left \lfloor \frac{M}{d} \right \rfloor}{|i-j|}}}}$$
（最后一步是将$dd‘$用$d$代换了，所以才会有$d‘|d$的限制。
后面那一坨只和$\frac{N}{d},\frac{M}{d}$有关系的式子的值是：设$A=\max(\frac{N}{d},\frac{M}{d}), B=\max(\frac{N}{d},\frac{M}{d})$，则$\sum_{i=1}^{A}{\sum_{j = 1}^{B}{|i-j|}}=\frac{(A-1)A(A+1)}{3}+\frac{AB(B-A)}{2}$，这里懒得推导，存在更优美的式子。
而前面那一坨，令$f(n)=\sum_{d|n}{d\mu(d)}$，由积性函数的性质可知，它也是积性函数，我们可以考虑筛出它的值。$f(1)=1$。
若$n$为质数，则$f(n)=1*\mu(1)+n*\mu(n)=1-n$。
若$n$的最小质因数$p$只出现了一次，则$f(n)=f(\frac{n}{p})+f(\frac{n}{p})*\mu(p)=f(\frac{n}{p})f(p)$。
若$n$的最小质因数$p$出现了两次以上，则$f(n)=f(\frac{n}{p})+f(\frac{n}{p})*0=f(\frac{n}{p})$。
对于询问，容易看出$(\left \lfloor \frac{N}{d} \right \rfloor,\left \lfloor \frac{M}{d} \right \rfloor)$的取值最多只有$2\left \lfloor \sqrt{N} \right \rfloor+2\left \lfloor \sqrt{M} \right \rfloor$种，且取值对应的$d$是一段连续的区间，可以将$f(n)$预处理出前缀和，$O(\left \lfloor \sqrt{N} \right \rfloor+\left \lfloor \sqrt{M} \right \rfloor)$回答询问，总时间复杂度$O(\max(N,M)+T(\left \lfloor \sqrt{N} \right \rfloor+\left \lfloor \sqrt{M} \right \rfloor))$。

代码：

#include <cstdio>
#include <algorithm>
using namespace std;
const int maxn = 2000001, mod1 = 1000000007, mod2 = 1000000009;
int t, n, m, p, q, tot, prime[maxn], sf[maxn][2], ans[2];
bool vis[maxn];
void inc(int &x, int y, int mod)
{
    x += y;
    if(x >= mod)
        x -= mod;
}
int f(int A, int B, int mod)
{
    if(!A || !B)
        return 0;
    if(A > B)
        swap(A, B);
    int cnt1 = (long long)(A - 1) * A * (A + 1) / 3 % mod;
    int cnt2 = (long long)A * B * (B - A) / 2 % mod;
    inc(cnt1, cnt2, mod);
    return cnt1;
}
int main()
{
    sf[1][0] = sf[1][1] = 1;
    for(int i = 2; i < maxn; ++i)
    {
        if(!vis[i])
        {
            prime[tot++] = i;
            sf[i][0] = mod1 - i + 1;
            sf[i][1] = mod2 - i + 1;
        }
        for(int j = 0; j < tot && (long long)i * prime[j] < maxn; ++j)
        {
            vis[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                sf[i * prime[j]][0] = sf[i][0];
                sf[i * prime[j]][1] = sf[i][1];
                break;
            }
            else
            {
                sf[i * prime[j]][0] = (long long)sf[i][0] * sf[prime[j]][0] % mod1;
                sf[i * prime[j]][1] = (long long)sf[i][1] * sf[prime[j]][1] % mod2;
            }
        }
    }
    for(int i = 2; i < maxn; ++i)
    {
        inc(sf[i][0], sf[i - 1][0], mod1);
        inc(sf[i][1], sf[i - 1][1], mod2);
    }
    scanf("%d", &t);
    while(t--)
    {
        scanf("%d%d", &n, &m);
        if(n > m)
            swap(n, m);
        ans[0] = ans[1] = 0;
        for(int i = 1, j; i <= n; i = j + 1)
        {
            j = min(n / (n / i), m / (m / i));
            inc(ans[0], (long long)(sf[j][0] - sf[i - 1][0] + mod1) * f(n / i, m / i, mod1) % mod1, mod1);
            inc(ans[1], (long long)(sf[j][1] - sf[i - 1][1] + mod2) * f(n / i, m / i, mod2) % mod2, mod2);
        }
        printf("%d %d\n", ans[0], ans[1]);
    }
    return 0;
}