SDOI 2013 项链

题面
首先，求出不同的珠子的数量。按照题意，珠子数量为：
${\sum }_{i=1}^a{\sum }_{j=i}^a{\sum }_{k=j}^a[gcd(i,j,k)=1]$
进行简单容斥，得到：
$\frac{1}{6}(2+{\sum }_{i=1}^a{\sum }_{j=1}^a(3[gcd(i,j)=1]+{\sum }_{k=1}^a[gcd(i,j,k)=1]))$
进行莫比乌斯反演，得到：
$\frac{1}{6}(2+{\sum }_{i=1}^a\mu (i)[\frac{a}{i}{]}^2([\frac{a}{i}]+3))$
线性筛预处理并整除分块即可。
设不同的珠子有 $m$ 种，设大小为 $d$ 不考虑旋转的项链数为 $f(n)$，考虑倒数第二个珠子是否与项链开头的珠子相同，则有递推式：
$f(1)=1,f(2)=m(m-1),f(n)=(m-2)f(n-1)+(m-1)f(n-2)$
求出 $f$ 的通项：
$f(n)=(m-1{)}^n+(-1{)}^{n}m$
根据Burnside引理，答案为：
$\frac{1}{n}{\sum }_{i=1}^{n}f(gcd(i,n))$
枚举最大公约数，得到答案为：
$\frac{1}{n}{\sum }_{d|n}f(\frac{n}{d})\phi (d)$
因此利用PollardRho算法分解 $n$ 的质因数，然后DFS枚举 $n$ 的所有约数，求出答案。
由于 $n$ 可能为模数的倍数，没有逆元，因此将任意整数记为 $kn+r$ 的形式，其中 $k$ 为整数，对模数取模。
时间复杂度 $O(t(a^{\frac{1}{2}}+n^{\frac{1}{4}}lo{g}_{2}n))$，空间复杂度 $O(a)$

#include<iostream>
#include<vector>
#include<algorithm>
#include<cmath>
#include<time.h>
using namespace std;
#define R register int
#define L long long
#define I inline
#define N 10000001
#define P 1000000007
int mu[N],prime[664579],b[6]={2,3,5,7,11,13},d[11];
bool vis[N];
L ps[11];
I L Add(L x,L y,L&p){
	L res=x;
	res+=y;
	return res>=p?res-p:res;
}
I L MulMod(L x,L y,const L p){
	L t=x*y-(L)((long double)x*y/p+0.5)*p;
	return t<0?t+p:t;
}
struct Number{
	int A;
	L B;
	I void Init(L x,L&n){
		A=x/n%P;
		B=x%n;
	}
};
I Number Addition(Number x,Number y,L&n){
	Number res;
	res.A=x.A+y.A;
	res.B=x.B+y.B;
	if(res.B>=n){
		res.B-=n;
		res.A++;
	}
	if(res.A>=P){
		res.A-=P;
	}
	return res;
}
I Number Sub(Number x,Number y,L&n){
	Number res;
	res.A=x.A-y.A;
	res.B=x.B-y.B;
	if(res.B<0){
		res.B+=n;
		res.A--;
	}
	if(res.A<0){
		res.A+=P;
	}
	return res;
}
I Number Mul(Number x,Number y,L&n){
	Number res;
	res.B=MulMod(x.B,y.B,n);
	res.A=(n%P*x.A%P*y.A+(L)y.B%P*x.A+x.B%P*y.A+(L)((long double)x.B*y.B/n))%P;
	return res;
}
I Number Pow(Number x,L y,L&n){
	Number res;
	res.Init(1,n);
	while(y!=0){
		if((y&1)==1){
			res=Mul(res,x,n);
		}
		y>>=1;
		x=Mul(x,x,n);
	}
	return res;
}
I Number GetF(L x,Number&m,L&n){
	Number res;
	res=Pow(m,x,n);
	if((x&1)==1){
		return Sub(res,m,n);
	}
	return Addition(res,m,n);
}
I Number Read(L C){
	int r,n;
	cin>>n;
	C*=6;
	Number res;
	res.Init(2,C);
	for(R i=1;i<=n;i=r+1){
		int t=n/i;
		r=n/t;
		Number tem,tem2;
		tem.Init((L)t*t,C);
		tem2.Init(t+3,C);
		tem=Mul(tem,tem2,C);
		if(mu[r]>mu[i-1]){
			tem2.Init(mu[r]-mu[i-1],C);
			res=Addition(res,Mul(tem,tem2,C),C);
		}else{
			tem2.Init(mu[i-1]-mu[r],C);
			res=Sub(res,Mul(tem,tem2,C),C);
		}
	}
	res.B/=6;
	return res;
}
I L PowMod(L x,L y,const L p){
	L res=1;
	while(y!=0){
		if((y&1)==1){
			res=MulMod(res,x,p);
		}
		y>>=1;
		x=MulMod(x,x,p);
	}
	return res;
}
I bool CheckPrime(const L n){
	if(n==2||n==3){
		return true;
	}
	if(n==1||n%6!=1&&n%6!=5){
		return false;
	}
	L p=n-1;
	int q=0;
	do{
		q++;
		p>>=1;
	}while((p&1)==0);
	for(R i=0;i!=6&&b[i]<n;i++){
		L t=PowMod(b[i],p,n);
		for(R j=0;j!=q&&t!=1;j++){
			L v=MulMod(t,t,n);
			if(v==1&&t!=n-1){
				return false;
			}
			t=v;
		}
		if(t!=1){
			return false;
		}
	}
	return true;
}
I L Gcd(L x,L y){
	return y==0?x:Gcd(y,x%y);
}
I L GetRand(){
	return (L)rand()<<45|(L)rand()<<30|rand()<<15|rand();
}
I L GoNext(L cur,const L c,L&n){
	return Add(MulMod(cur,cur,n),c,n);
}
I L Pollard(L n){
	do{
		L c=GetRand()%(n-1)+1,t,r;
		t=c;
		r=GoNext(t,c,n);
		while(t!=r){
			L d=Gcd(t-r,n);
			if(d<0){
				d=-d;
			}
			if(d!=1){
				return d;
			}
			t=GoNext(t,c,n);
			r=GoNext(GoNext(r,c,n),c,n);
		}
	}while(true);
}
I void GetPrime(L n,vector<L>&D){
	if(CheckPrime(n)==true){
		return D.push_back(n);
	}
	int t=sqrt(n);
	if((L)t*t==n){
		return GetPrime(t,D);
	}
	L d=Pollard(n);
	GetPrime(d,D);
	GetPrime(n/d,D);
}
I int GetD(L n,L p){
	int res=0;
	do{
		res++;
		n/=p;
	}while(n%p==0);
	return res;
}
I void DFS(int t,L s,L phi,Number&m,L&n,Number&ans){
	if(t==-1){
		Number tem;
		tem.Init(phi,n);
		ans=Addition(ans,Mul(tem,GetF(n/s,m,n),n),n);
		return;
	}
	DFS(t-1,s,phi,m,n,ans);
	phi*=ps[t]-1;
	for(R i=0;i!=d[t];i++){
		s*=ps[t];
		DFS(t-1,s,phi,m,n,ans);
		phi*=ps[t];
	}
}
I void Solve(){
	L n;
	cin>>n;
	int t=0;
	Number m=Read(n),e;
	e.Init(1,n);
	m=Sub(m,e,n);
	if(n==1){
		cout<<(n%P*m.A+m.B)%P<<endl;
		return;
	}
	vector<L>D;
	GetPrime(n,D);
	sort(D.begin(),D.end());
	ps[0]=D[0];
	d[0]=GetD(n,ps[0]);
	for(vector<L>::iterator T=D.begin()+1;T!=D.end();T++){
		if(*T!=*(T-1)){
			t++;
			ps[t]=*T;
			d[t]=GetD(n,*T);
		}
	}
	Number ans;
	ans.Init(0,n);
	DFS(t,1,1,m,n,ans);
	cout<<ans.A<<endl;
}
int main(){
	srand((int)time(0));
	mu[1]=1;
	int t=0;
	for(R i=2;i!=N;i++){
		if(vis[i]==false){
			prime[t]=i;
			t++;
			mu[i]=-1;
		}
		for(R j=0;prime[j]*i<N;j++){
			vis[i*prime[j]]=true;
			if(i%prime[j]==0){
				mu[i*prime[j]]=0;
				break;
			}
			mu[i*prime[j]]=-mu[i];
		}
		mu[i]+=mu[i-1];
	}
	cin>>t;
	for(R i=0;i!=t;i++){
		Solve();
	}
	return 0;
}