[HDU 5780] gcd (公式证明)

做出这题你需要推出一个重要的式子：$gcd(x^a-1,x^b-1)=x^{gcd(a,b)}-1$
我这证明可能不算严谨吧。。。。
反正OI不需要证明，只需要感性理解。然而我个人觉得感性理解反而比证明重要啊，证明不就是几个式子套来套去，过几天就忘光了。
不妨设$a>b$，利用根相减损术，可以得到 $$gcd(x^a-1,x^b-1)=gcd(x^a-x^b,x^b-1)$$
化简得 $$gcd(x^a-1,x^b-1)=gcd(x^b(x^{a-b}-1),x^b-1)$$
因为 $$gcd(x^b,x^b-1)=1$$
所以$$gcd(x^a-1,x^b-1)=gcd(x^{a-b}-1,x^b-1)$$
至此，我们可以发现 $a$ 和 $b$ 也恰好进行了一次根相减损术的过程。
根相减损术结束条件是其中一边为 $0$，另一边就是 $gcd$ 。
不妨设最后$a=0$，则 $b$ 就为原来 $a,b$ 的 $gcd$
带入原式，得$$gcd(x^a-1,x^b-1)=gcd(x^0-1,x^{gcd(a,b)}-1)$$
就是$$gcd(x^a-1,x^b-1)=x^{gcd(a,b)}-1$$
于是就证完了，题目就可以转化为 $[1,n]$ 中每对数的 $gcd$ 对答案的贡献
首先可以想到的思路是枚举 $d=gcd(a,b)$，有序对$(a,b)$对数就是${\sum }_{i=1}^{n/d}\phi (i)$（这就相当于找出所有互质的数再同时乘上 $gcd$），可以把 $\phi (i)$ 前缀和预处理出来。题目求的是无序对，就乘2再把$a,b$相同的情况减掉（因为a,b相同的情况算了两次）。
我就只想到这儿了，但到这里还是$O(300\ast n)$，虽然看起来还挺优但肯定过不了。
其实只差最后一步了，看到 $n/d$ 就弄个分块嘛，然后还有个等差数列求和，可是偏偏这步没想出来QAQ太菜了QAQ
还有需要注意的一点是特判$x=1$的情况，还有$x^{gcd(a,b)}-1$里的那个1别忘了减。

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1000001;
const int p=1e9+7;
int T,x,n;
int phi[N];
bool b[N];

void read(int &x){
	char ch=getchar();x=0;
	for(;ch<'0'||ch>'9';ch=getchar());
	for(;ch>='0'&&ch<='9';ch=getchar()) x=(x<<3)+(x<<1)+ch-'0';
}

ll power(int a,int n){
	ll ans=1;
	for(ll sum=a;n;sum=sum*sum%p,n>>=1) if (n&1) ans=ans*sum%p;
	return ans;
}

void Add(ll &x,int y){
	x+=y;
	while(x<0) x+=p;
	while(x>=p) x-=p;
}

void Add(int &x,int y){
	x+=y;
	while(x<0) x+=p;
	while(x>=p) x-=p;
}

void init(){
	for(int i=1;i<N;i++) phi[i]=i;
	for(int i=2;i<N;i++)
	 if (!b[i]){
	 	phi[i]=i-1;
	 	for(int j=2;i*j<N;j++)
	 	 b[i*j]=1,phi[i*j]=phi[i*j]/i*(i-1);
	 }
	for(int i=2;i<N;i++) Add(phi[i],phi[i-1]);
}

int main(){
	read(T);
	init();
	for(int o=1;o<=T;o++){
		read(x);read(n);
		if (x==1){ printf("0\n");continue; }
		ll ans=0,inv=power(x-1,p-2);
		for(int i=1,j=1;i<=n;i=j+1){
			j=n/(n/i);
			int w=n/i;
			ll v=power(x,i)*(power(x,j-i+1)-1)%p*inv%p*phi[w]%p*2%p;
			Add(ans,(int)v);
			//cout<<i<<' '<<j<<' '<<power(x,i)<<' '<<(power(x,j-i+1)-1)*inv%p<<endl;
		}
		Add(ans,-1LL*n*n%p);
		ans=(ans-(power(x,n)-1)*x%p*inv%p)%p;
		printf("%lld\n",(ans+p)%p);
	}
	return 0;
}