题意
求 $f(n,a,b)=sum_{i=1}^n sum_{j=1}^i gcd(i^a-j^a,i^b-j^b)[gcd(i,j)=1]\%(10^9+7)$，$1 le n,a,b le 10^9$，共有 $T$ 组测试，其中只有10组的 $n$ 大于 $10^6$.
分析
首先，当 $i, j$互质，$a, b$互质时，有 $gcd(i^a-j^a,i^b-j^b)=i-j$（证明见 链接），也可以打表猜一猜嘛。
可以推出：$$sum_{d=1}^{N}mu(d)cdot dsum_{i=1}^{lfloorfrac{N}{d} floor}sum_{j=1}^{i}(i-j)$$
单独考虑后半部分，$sum_{i=1}^{k}sum_{j=1}^{i}(i-j)=frac{k^3-k}{6}$.
然后，只剩下左边的 $mu(d)cdot d$,
将其与恒等函数 $Id(n) = n$ 狄利克雷卷积后得
$$egin{align*}
(mu(d)cdot d)*Id(d)
& = sum_{d|n}(mu(d)cdot d)cdot Id(frac{n}{d})\
& = sum_{d|n}mu(d) = [n=1]
end{align*}$$
接下来套杜教筛得公式
$$egin{align*}
S(n)
& = sumlimits_{i=1}^n [i=1]-sumlimits_{i=2}^nicdot S(lfloordfrac{n}{i}
floor)\
& = 1-sumlimits_{i=2}^nicdot S(lfloordfrac{n}{i}
floor)
end{align*}$$
#include<bits/stdc++.h> using namespace std; typedef long long ll; const int maxn = 6e6 + 10; const ll mod = 1e9+7; const ll inv6 = 166666668; int sum[maxn], mu[maxn], pri[maxn], pn; bool vis[maxn]; map<int, int>mp_sum; int n, a, b; ll s2(ll i, ll j) { return (i+j) * (j-i+1) / 2 % mod; } ll s3(ll k) { return (k*k%mod - 1) * k % mod * inv6 % mod; } ll S(ll x) { if(x < maxn) return sum[x]; if(mp_sum[x]) return mp_sum[x]; ll ret = 1LL; for(int i = 2, j;i <= x;i = j+1) { j = x / (x / i); ret = (ret - s2(i, j) * (S(x/i))%mod) % mod; } return mp_sum[x] = (ret + mod) % mod; } void pre() { mu[1] = 1; for(int i = 2;i < maxn;i++) { if(!vis[i]) { pri[++pn] = i; mu[i] = -1; } for(int j = 1;j <= pn && i * pri[j] < maxn; j++) { vis[i * pri[j]] = true; if(i % pri[j]) mu[i * pri[j]] = -mu[i]; else { mu[i * pri[j]] = 0; break; } } } for(int i = 1;i < maxn;i++) sum[i] = (sum[i-1] + i * mu[i]) % mod; } int main() { pre(); int T; scanf("%d", &T); while(T--) { scanf("%d%d%d", &n, &a, &b); ll ans = 0; for(ll l = 1,r; l <= n;l = r+1) { r = n / (n / l); ans = (ans + (S(r) - S(l-1)) * s3(n/l)) % mod; } printf("%lld ", (ans+mod)%mod); } return 0; }
最开始开的 MAXN=2e6，会TLE；原博客开的6e6，又MLE，将long long 数组改成 int 才行。
其实标答是推成 $displaystyle ans = frac{sum _{i=1}^n ivarphi (i) - 1}{2}$，少了一次整除分块。
但是，通过这种解法，让我深刻认识了杜教筛的时空矛盾该怎么平衡。