HDU 6051 If the starlight never fade（原根+推式子）

题目大意：

设\(f(i)\)为使\((x+y)^i \equiv x^i (mod\ p)\)成立的(x,y)的对数。其中\(1 \leq x \leq p-1 , 1\leq y\leq m\),m，p给定且p是一个质数。求\(\sum_{i=1}^{p-1}i*f(i)\),p<=1e9+7,m<=p-1

思路

我们考虑用原根去代换x，y。
设g为p的一个原根，\(g^a\equiv x(mod \ p),g^b \equiv y(mod \ p)\)。
然后我们用\(g\)去代换\(x,y\)。\((g^a+g^b)^i\equiv g^{a*i}(mod\ p)\)然后我们在式子两边同时除一个\(g^{a*i}\)。得到\((1+g^{b-a})^i\equiv 1(mod\ p)\)
设\(1+g^{b-a}\equiv g^k(mod\ p)\)(因为原根的性质所以我们一定可以找到这样的k)。
此时原式为\(g^{k*i} \equiv 1(mod\ p)\)由费马小定理可以到到\(p-1 \mid k*i\)。
这就要求k为\(\frac{p-1}{gcd(p-1,i)}\)的倍数(即至少包含p-1所特有的因子)。由于\(0<k<p-1\),为什么k不能取0呢？因为$g^{b-a}不会为0，所以\(1+g^{b-1}> 1\)。
所以可以得到这样的k的数量是\(\frac{p-1}{\frac{p-1}{gcd(p-1,i)}}-1=gcd(p-1,i)-1\)，这里因为k不能取0和p-1，所以要减1。
又因为\(1+g^{b-a} \equiv g^k(mod\ p)\)得\(g^b \equiv(g^k-1)*g^a (mod\ p)\)，\(y\equiv x*(g^k-1)(mod\ p)\)。每有一个k，y对应一个x。所以对于一个y有\(gcd(p-1,i)-1\)个x对应。
所以可以得出\(f(i)=m*(gcd(p-1,i)-1)\)。
\(\sum_{i=1}^{p-1}i*f(i)\)
\(=m\sum_{i=1}^{p-1}i*gcd(p-1,i)-m*\frac{(p-1)*p}{2}\)
重点是如何求\(\sum_{i=1}^{p-1}i*gcd(p-1,i)\)。
\(\sum_{i=1}^{p-1}i*gcd(p-1,i)\)
显然d是p-1的约数
\(=\sum_{d\mid p-1}d\sum_{i=1}^{p-1}i*[gcd(p-1,i)==d]\)
\(=\sum_{d\mid p-1}d^2\sum_{i=1}^{\frac{p-1}{d}}i*[gcd(\frac{p-1}{d},i)==1]\)
然后有一个神奇的变换。
\(\sum_{i=1}^{t}i*[gcd(t,i)==1]=\frac{t*\varphi(t)+[t==1]}{2}\)
为什么？
令\(t=\frac{p-1}{d}\)就是\(\sum_{i=1}^{t}i*[gcd(t,i)==1]\)
其实就是求1到t中与t互质的数的和。
由更相减损术得若\(gcd(n,i)=1\)则\(gcd(n,n-i)=1\)
所以一个与\(t\)互质的数\(x\)，\(t-x\)也与\(t\)互质。
所以与t互质的数成对出现，设这一对数为a,b，有\(a+b=t\)。
所以\(\sum_{i=1}^{t}i*[gcd(t,i)==1]=\frac{t*\varphi(t)+[t==1]}{2}\)
故原式可以化为
\(\sum_{d\mid p-1}d^2\sum_{i=1}^{\frac{p-1}{d}}i*[gcd(\frac{p-1}{d},i)==1]\)
\(=\sum_{d\mid p-1}d^2*\frac{\frac{p-1}{d}*\varphi(\frac{p-1}{d})+[\frac{p-1}{d}==1]}{2}\)
带回去求解即可。
复杂度？\(O(能过)\)，因为求约数可以先扫出质数来优化，所以\(\sqrt{\frac{p}{lnp}}*\sqrt{p}=\sqrt{\frac{p^2}{lnp}}\)差不多两亿，不是很大。而且求欧拉函数是求p-1的约数的约数，一定比\(\sqrt{p}\)小，所以能过？

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
#define int long long
const int N=1001000;
const int mod=1e9+7;
bool book[N];
int prime[N],cnt,inv,T;
int read(){
	int sum=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9'){sum=sum*10+ch-'0';ch=getchar();}
	return sum*f;
}
void pre_work(){
	for(int i=2;i<=1000000;i++){
		if(book[i]==0)prime[++cnt]=i;
		for(int j=1;j<=cnt&&prime[j]*i<=1000000;j++){
			book[prime[j]*i]=1;
			if(i%prime[j]==0)break;
		}
	}
}
int ksm(int x,int b){
	int tmp=1;
	while(b){
		if(b&1)tmp=tmp*x%mod;
		b>>=1;
		x=x*x%mod;
	}
	return tmp;
}
int phi(int x){
	int tmp=x;
	int ans=x;
	for(int i=1;i<=cnt&&prime[i]*prime[i]<=x;i++){
		if(tmp%prime[i]==0){
			ans=ans/prime[i]*(prime[i]-1);
			while(tmp%prime[i]==0)tmp/=prime[i];
		}
	}
	if(tmp>1)ans=ans/tmp*(tmp-1);
	return ans;
}
int work(int x){
	int tmp=0;
	for(int i=1;i*i<=x;i++){
		if(x%i==0){
			int a=i*i%mod;
			int b=x/i;
			int c=phi(b);
			int d=(b==1);
			tmp=(tmp+a*(b*c%mod+d)%mod*inv%mod)%mod;
			if(x/i>i){
				int hh=x/i;
				int a=hh*hh%mod;
				int b=x/hh;
				int c=phi(b);
				int d=(b==1);
				tmp=(tmp+a*(b*c%mod+d)%mod*inv%mod)%mod;
			}
		}
	}
	return tmp;
}
signed main(){
	T=read();
	pre_work();
	inv=ksm(2,mod-2);
	int now=0;
	while(T--){
		now++;
		int m=read(),p=read();
		int tmp=((m*work(p-1)%mod-m*(p-1ll)%mod*p%mod*inv%mod)%mod+mod)%mod;
		printf("Case #%lld: %lld\n",now,tmp);
	}
	return 0;
}