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