BZOJ 4018 小Q的幻想之乡

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

题意：
T组询问，每组询问给定N和M，求∑Ni=1∑Mj=1|i−j|gcd(i,j)分别模109+7和109+9的值。
T≤1000,N,M≤2∗106。

题解：
考虑将∑Ni=1∑Mj=1|i−j|gcd(i,j)化简。
首先引入一个记号，令[x]表示：若x为真，则[x]为1；否则[x]为0，其中x是一个布尔表达式。

∑i=1N∑j=1M|i−j|gcd(i,j)=∑i=1N∑j=1M∑d|i−j|gcd(i,j)[d=gcd(i,j)]=∑d∑i=1⌊Nd⌋∑j=1⌊Md⌋|i−j|[gcd(i,j)=1]=∑d∑i=1⌊Nd⌋∑j=1⌊Md⌋|i−j|∑d′|gcd(i,j)μ(d′)=∑d∑d′d′μ(d′)∑i=1⌊Ndd′⌋∑j=1⌊Mdd′⌋|i−j|=∑d∑d′|dd′μ(d′)∑i=1⌊Nd⌋∑j=1⌊Md⌋|i−j|

（最后一步是将dd′用d代换了，所以才会有d′|d的限制。
后面那一坨只和Nd,Md有关系的式子的值是：设A=min(Nd,Md),B=max(Nd,Md)，则∑Ai=1∑Bj=1|i−j|=(A−1)A(A+1)3+AB(B−A)2，这里懒得推导，存在更优美的式子。
而前面那一坨，令f(n)=∑d|ndμ(d)，由积性函数的性质可知，它也是积性函数，我们可以考虑筛出它的值。f(1)=1。
若n为质数，则f(n)=1∗μ(1)+n∗μ(n)=1−n。
若n的最小质因数p只出现了一次，则f(n)=f(np)+f(np)∗μ(p)=f(np)f(p)。
若n的最小质因数p出现了两次以上，则f(n)=f(np)+f(np)∗0=f(np)。
对于询问，容易看出(⌊Nd⌋,⌊Md⌋)的取值最多只有2⌊N−−√⌋+2⌊M−−√⌋种，且取值对应的d是一段连续的区间，可以将f(n)预处理出前缀和，O(⌊N−−√⌋+⌊M−−√⌋)回答询问，总时间复杂度O(max(N,M)+T(⌊N−−√⌋+⌊M−−√⌋))。

代码：

#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;}