HDU 6706

2019CCPC网络赛1005
问题求解：
$\sum _{i=1}^n\sum _{j=1}^{i}gcd(i^a-j^a,i^b-j^b)[gcd(i,j)=1]\%(mod),mod=1e9+7,a和b互素，且a\ge b$
引理1：$x^n-y^n=(x-y)(x^{n-1}+x^{n-2}y+\dots +y^{n-1})$
引理2：当$i,j$互素时，$gcd(i^a-j^a,i^b-j^b)=i^{gcd(a,b)}-j^{gcd(a,b)}$
证明：记$a=k\ast b+r,$
$i^a-j^a=i^{k\ast b+r}-i^{k\ast b}\ast j^r+i^{k\ast b}\ast j^r-j^{k\ast b+r}$
$=i^{k\ast b}\ast (i^r-j^r)+j^r\ast (i^{k\ast b}-j^{k\ast b})$
$=j^r(i^b-j^b)(i^{b\ast (k-1)}+i^{b\ast (k-2)}j+\dots +j^{b\ast (k-1)})+i^{k\ast b}(i^r-j^r)$
$gcd(i^a-j^a,i^b-j^b)$（引理1）
$=gcd(i^{k\ast b}(i^r-j^r),i^b-j^b)$
因为$i^{kb}-j^{kb}=(i^b-j^b)(i^{b\ast (k-1)}+i^{b\ast (k-2)}j+\dots +j^{b\ast (k-1)})$
所以记$d=gcd(i^{kb},i^b-j^b)=gcd(j^{kb},i^b-j^b)$
则$d|i^{kb},d|j^{kb},d|i^b-j^b$,因为$i,j$互素，所以$d=1$
$gcd(i^a-j^a,i^b-j^b)=gcd(i^{k\ast b}(i^r-j^r),i^b-j^b)$
$=gcd(i^r-j^r,i^b-j^b)=i^{gcd(a,b)}-j^{gcd(a,b)}$(类似辗转相除法)
证毕。
引理3：$\sum _{i=1}^n\sum _{j=1}^{i}j[gcd(i,j)=1]=\sum _{i=1}^n(\frac{i\phi (i)}{2}+[i=1])$
证明：$gcd(a,x)=gcd(x-a,x),a\le x$
若$a$和$x$互素，则$x-a$和$x$也互素，所以
$\sum _{i=1}^n\sum _{j=1}^{i}j[gcd(i,j)=1]=\frac{\sum _{i=1}^n\sum _{j=1}^{i}i[gcd(i,j)=1]}{2}+[i=1]$
$=\sum _{i=2}^n(\frac{i\phi (i)}{2})+1$
证毕。
$\sum _{i=1}^n\sum _{j=1}^{i}gcd(i^a-j^a,i^b-j^b)[gcd(i,j)=1]$
$=\sum _{i=1}^n\sum _{j=1}^i(i^{gcd(a,b)}-j^{gcd(a,b)})[gcd(i,j)=1]$（引理2）
$=\sum _{i=1}^n\sum _{j=1}^i(i-j)[gcd(i,j)=1]$,($a,b$互素)
$=\sum _{i=1}^n\sum _{j=1}^{i}i[gcd(i,j)=1]-\sum _{i=1}^n\sum _{j=1}^{i}j[gcd(i,j)=1]$（引理3）
$=\sum _{i=2}^n(\frac{i\phi (i)}{2})$
$=\frac{\sum _{i=1}^{n}i\phi (i)-1}{2}$
现在问题转换为了如何求解$\sum _{i=1}^{n}i\phi (i)$,$n\le 1e9$
构造杜教筛：
$f=id\cdot \phi ,g=id,h=f\ast g=i^2$
$\sum _{i=1}^{n}h(i)=\sum _{i=1}^n\sum _{d|i}f(d)g(\frac{i}{d})$
$=\sum _{i=1}^n\sum _{d|i}d\phi (d)\cdot \frac{i}{d}$
$=\sum _{i=1}^{n}i\sum _{d|i}\phi (d)=\sum _{i=1}^n{i}^2$
利用杜教筛公式：$g(1)S(n)=\sum _{i=1}^{n}h(i)-\sum _{d=2}^{n}g(d)S(\lfloor \frac{n}{d}\rfloor )$
后面就是杜教筛了。
AC代码：

#include<bits\stdc++.h>
#define ll long long
using namespace std;
const int mod = 1e9+7;
const int maxn = 1e6+10;
ll p[maxn];
ll sum[maxn];
ll inv2,inv6;
map<int,ll>mp;
ll qpow(ll a,ll b){
	ll res=1;
	while(b!=0){
		if(b&1) res=res*a%mod;
		a=a*a%mod;
		b>>=1;
	}
	return res;
}
ll f_g_sum(int n){
	return (ll)n%mod*(n+1)%mod*(2*n+1)%mod*inv6%mod;
}
ll g_sum(int n){
	return (ll)n%mod*(n+1)%mod*inv2%mod;
}
ll GetSum(int n){
	if(n<=1000000)return sum[n];
	if(mp[n]!=0)return mp[n];
	ll ans = f_g_sum(n)%mod; 
	for(ll l = 2, r; l <= n; l = r + 1) { 
		r = (n / (n / l)); 
		ans = (ans - (g_sum(r) - g_sum(l - 1)) * GetSum(n / l) + mod)%mod;
	}
	mp[n]=ans;
	return ans; 
}
void phi_table(int n, ll* phi) {
  for (int i = 2; i <= n; i++) phi[i] = 0;
  phi[1] = 1;
  for (int i = 2; i <= n; i++)
    if (!phi[i])
      for (int j = i; j <= n; j += i) {
        if (!phi[j]) phi[j] = j;
        phi[j] = phi[j] / i * (i - 1);
      }
}
int main(){
	int T,n,a,b;
	scanf("%d",&T);
	inv2=qpow(2,mod-2);
	inv6=qpow(6,mod-2);
	phi_table(1000000,p);
	sum[0]=0;
	for(int i=1;i<=1000000;i++){
		sum[i]=(sum[i-1]+i*p[i])%mod;
	}
	while(T--){
		scanf("%d%d%d",&n,&a,&b);
		ll ans;
		ans=(GetSum(n)-1+mod)%mod*inv2%mod;
		printf("%lld\n",ans);
	}
	return 0;
}