打表观察得到，gcd(i,j)==1时，gcd(i^a−j^a,i^b−j^b)的值为i - j。所以，你发现这个题跟ab就没关系了...
变成去求∑∑(i-j)[gcd(i,j) == 1]了。有一个显然的结论，gcd(i,j) == gcd(i-j,i)。
设k为i-j，则变成
∑(i 1->n)∑k(1->i-1)[gcd(i,k) == 1]。
又因为i和i本身，一定不互质，所以变成
∑(i 1->n)∑k(1->i)[gcd(i,k) == 1]。然后有个结论，1…N中与N互质的数的和，ans=N*phi(N)/2。
然后就变成了求∑phi(i)*i，求前缀和想到了杜教筛。卷积一下函数g(x) = x，就行了。
1 #include <cstdio> 2 #include <map> 3 #include <cmath> 4 using namespace std; 5 typedef long long ll; 6 const int MAXN = 1000100,mo = 1e9 + 7,inv2 = (mo + 1) / 2,inv6 = (mo + 1) / 6; 7 8 int T,n,a,b,maxn,phi[MAXN],sum[MAXN],pri[MAXN]; 9 bool vis[MAXN]; 10 map<int,int> f; 11 int solve(int n) 12 { 13 if (n <= maxn) 14 return sum[n]; 15 if (f.count(n)) 16 return f[n]; 17 int ans = (ll)n * (n + 1) % mo * (n * 2 + 1) % mo * inv6 % mo; 18 for (int l = 2,r;l <= n;l = r + 1) 19 { 20 r = n / (n / l); 21 int tp = (ll)(l + r) * (r - l + 1) % mo * inv2 % mo; 22 ans -= (ll)tp * solve(n / l) % mo; 23 if (ans < 0) 24 ans += mo; 25 } 26 return f[n] = ans; 27 } 28 void init() 29 { 30 maxn = 1000000; 31 phi[1] = 1; 32 int tot = 0; 33 for (int i = 2;i <= maxn;i++) 34 { 35 if (vis[i] == false) 36 { 37 pri[++tot] = i; 38 phi[i] = i - 1; 39 } 40 for (int j = 1;j <= tot && i * pri[j] <= maxn;j++) 41 { 42 vis[i * pri[j]] = true; 43 if (i % pri[j] != 0) 44 phi[i * pri[j]] = phi[i] * (pri[j] - 1); 45 else 46 { 47 phi[i * pri[j]] = phi[i] * pri[j]; 48 break; 49 } 50 } 51 } 52 for (int i = 1;i <= maxn;i++) 53 sum[i] = (sum[i - 1] + (ll)phi[i] * i % mo) % mo; 54 } 55 56 int main() 57 { 58 init(); 59 for (scanf("%d",&T);T != 0;T--) 60 { 61 scanf("%d%d%d",&n,&a,&b); 62 printf("%d ",(ll)(solve(n) - 1 + mo) % mo * inv2 % mo); 63 } 64 return 0; 65 }