【BZOJ3202】项链（莫比乌斯反演，Burnside引理）
题面
题解
首先读完题目，很明显的感觉就是，分成了两个部分计算。
首先计算本质不同的珠子个数，再计算本质不同的项链个数。
前面一个部分和(gcd)相关，一种莫比乌斯反演的感觉。
后面一个部分出现了旋转操作，要求本质不同的方案数，不难想到Burnside引理。
首先先考虑怎么求本质不同的珠子个数。
我们直接考虑无序的三元组((x,y,z))，满足(x,y,zle a,gcd(x,y,z)=1)
容斥考虑最终答案，就是(frac{1}{6}(S3-S2*3+S1*2))，其中(S3)是无序的三元组个数，(S2)是满足互质的二元组个数，而(S1)则是三个数都相同并且(gcd=1)的三元组个数。
(S1)显然只有一个，即可以构成唯一合法三元组((1,1,1))
考虑怎么计算(S2)和(S3)。
看起来(s2)很熟悉，也就是(sum_{i=1}^asum_{j=1}^a[gcd(i,j)=1])
然后直接莫比乌斯反演。
设(F(x)=sum_{i=1}^asum_{j=1}^a[gcd(i,j)=x])
(G(x)=sum_{x|d}F(d)=sum_{i=1}^asum_{j=1}^a[x|gcd(i,j)]=[frac{a}{x}]^2)
然后又因为(F(x)=sum_{x|d}G(d)*mu(frac{d}{x}))
要求的是(F(1))，即(F(1)=sum_{i=1}^aG(i)*mu(i)=sum_{i=1}^a[frac{a}{i}]^2mu(i))。
然后是(S3)，即三元组个数，和上面一样的推法，可以得到贡献是(sum_{i=1}^a[frac{a}{i}]^3mu(i))
好了，那么计算本质不同的珠子的个数就算完了。
接下来考虑如何计算本质不同的项链的个数。
因为只有旋转计算本质不同，那么置换一共(n)个，是顺时针旋转(i)位。
那么答案就是所有置换的不动点个数除以(n)。那么不难往Polya定理方面靠。
只需要知道每个循环的大小即可。
对于顺时针旋转(i)位的循环长度显然就是(frac{n}{gcd(n,i)})，循环个数(gcd(n,i))也很显然。
因为是不动点的个数，因此每一个循环内的所有珠子的颜色必须相同。同时，题目限制相邻两个位置的珠子颜色不能相同。因此，这里可以等价的看做有(gcd(n,i))个珠子，收尾相连，要求相邻的珠子不同色的方案数，假装这个的方案数是(f(x)),表示(x)个珠子收尾相连时候的方案数。那么现在要求的答案就是(sum_{i=1}^nf(gcd(i,n)))。(n)太大了，考虑优化。
我们枚举(gcd),即(sum_{d=1}^nsum_{i=1}^n[gcd(i,n)=d]f(d))，那么显然(d)是(n)的因数。所以可以接着写成(sum_{d|n}sum_{i=1}^n[gcd(i,n)=d]f(d)),进一步推，也就是(sum_{d|n}sum_{i=1,d|i}^n[gcd(frac{i}{d},frac{n}{d})=1]f(d))，即需要知道所有和(frac{n}{d})互质的(frac{i}{d})，考虑(frac{i}{d})的最大取值只有(frac{n}{d})，所以这个值显然就是(varphi(frac{n}{d}))。因此，总的贡献就是(sum_{d|n}f(d)varphi(frac{n}{d}))。
现在来考虑怎么计算(f(x)),假设一共有(m)中不同的珠子
首先我们可以在两个珠子之间插入一个不同的珠子，即(f(x-1)*(m-2))，亦或者强制让首位两个珠子颜色相同，然后插入一个不同颜色的珠子将他们隔开，这个操作等价于我现有(f(x-2))个珠子构成了一个环，然后插入一个和首位相同的珠子，再插入一个颜色不同的珠子隔开，即(f(x-2)*(m-1))，那么转移就是(f(x)=f(x-1)*(m-2)+f(x-2)*(m-1))。
显然可以矩阵乘法直接算。然而这种东西一般也可以推递推式的。
我们尽可能让他往一个等比的方向上靠，
那么我们最好能把式子化简为这种形式：(f(x)+af(x-1)=bf(x-1)+abf(x-2))
那么列出方程组：
所以有(b=a+m-2)，所以有(a(a+m-2)-(m-1)=0),也就是(a^2+(m-2)a-(m-1)=0)
因式分解一下就是((a+(m-1))(a-1)=0)
解出来(a=1)或者(a=1-m)，看着(a=1)好做些，
那么式子写成(f(x)+f(x-1)=(m-1)*(f(x-1)+f(x-2)))
所以得到(f(x)+f(x-1)=(m-1)^{x-2}(f(1)+f(2))),
显然有(f(1)=0,f(2)=m*(m-1))，为了方便，设(S(x)=f(x)+f(x-1))。
递推式换成(S(x)=(m-1)S(x-1)=(m-1)^{x-2}S(2)=(m-1)^{x-1}m)
考虑一下(f(x))怎么计算，(f(x)=S(x)-f(x-1)=m(m-1)^{x-1}-f(x-1))
一般这样子都可以构建一个递推式。
所以，我们可以构建出这样一个等式：
在考虑一下边界情况，可以得到(f(x)-(m-1)^x=(-1)^x(m-1))
即(f(x)=(-1)^x(m-1)+(m-1)^x),这样子就可以快速幂计算(f(x))了。
然而发现一些很不舒服的东西，因为(n)太大，导致(n)可能是模数的倍数，那就不能够直接除(n)了。
怎么办呢？那么首先我们直接模(mod^2)，这样子最终的结果会是(mod)的倍数，所以可以先除掉(mod)然后再考虑逆元就行了。
离线一下就跑得飞快了，目前是洛谷rk1，bzoj rk3
#include<iostream>
#include<cstdio>
#include<cmath>
using namespace std;
#define ll long long
const ll mod=1000000007;
const ll MOD=mod*mod;
const ll inv6 = 833333345000000041ll;
#define MAX 10000001
inline ll read()
{
ll x=0;bool t=false;char ch=getchar();
while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
if(ch=='-')t=true,ch=getchar();
while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
return t?-x:x;
}
void add(ll &x,ll y){x+=y;if(x>=MOD)x-=MOD;}
ll Multi(ll x,ll y){return (x*y-(ll)(((long double)x*y+0.5)/(long double)MOD)*MOD+MOD)%MOD;}
ll fpow(ll a,ll b)
{
ll s=1;
while(b){if(b&1)s=Multi(s,a);a=Multi(a,a);b>>=1;}
return s;
}
ll qpow(ll a,ll b)
{
ll s=1;
while(b){if(b&1)s=s*a%mod;a=a*a%mod;b>>=1;}
return s;
}
ll sqr(ll x){return Multi(x,x);}
ll cub(ll x){return Multi(sqr(x),x);}
bool zs[MAX];
int pri[MAX/10],tot,mu[MAX];
ll n,m,ans;int a,cnt;
ll nn[50],aa[50];
ll p[MAX/10],q[MAX/10];
void pre(ll N)
{
mu[1]=1;
for(int i=2;i<=N;++i)
{
if(!zs[i])pri[++tot]=i,mu[i]=-1;;
for(int j=1;j<=tot&&i*pri[j]<=N;++j)
{
zs[i*pri[j]]=true;
if(i%pri[j])mu[i*pri[j]]=-mu[i];
else break;
}
}
for(int i=1;i<=N;++i)mu[i]+=mu[i-1];
}
ll Calc(ll n)
{
ll S2=0,S3=0;
for(ll i=1,j;i<=n;i=j+1)
{
j=n/(n/i);
add(S3,Multi(cub(n/i),(mu[j]-mu[i-1]+MOD)%MOD));
add(S2,Multi(sqr(n/i),(mu[j]-mu[i-1]+MOD)%MOD));
}
add(S3,Multi(S2,3)),add(S3,2);
return Multi(S3,inv6);
}
ll F(ll n)
{
ll ret=fpow(m-1,n);
if(n&1)add(ret,(MOD-(m-1))%MOD);
else add(ret,m-1);
return ret;
}
void dfs(int i,ll d,ll phi)
{
if(i==cnt+1){add(ans,Multi(phi,F(n/d)));return;}
dfs(i+1,d,phi);
d*=p[i];phi*=p[i]-1;dfs(i+1,d,phi);
for(int x=2;x<=q[i];++x)
d*=p[i],phi*=p[i],dfs(i+1,d,phi);
}
int main()
{
int T=read();ll mx=0;
for(int i=1;i<=T;++i)nn[i]=read(),aa[i]=read();
for(int i=1;i<=T;++i)mx=max(mx,max((ll)sqrt(nn[i]),aa[i]));
pre(mx);
for(int TT=1;TT<=T;++TT)
{
n=nn[TT];a=aa[TT];
m=Calc(a);cnt=0;ll x=n;
for(int i=1;i<=tot&&1ll*pri[i]*pri[i]<=x;++i)
if(x%pri[i]==0)
{
p[++cnt]=pri[i];q[cnt]=0;
while(x%pri[i]==0)++q[cnt],x/=pri[i];
}
if(x>1)p[++cnt]=x,q[cnt]=1;
ans=0;dfs(1,1,1);
if(n%mod==0)ans=(ll)(ans/mod)*qpow(n/mod,mod-2)%mod;
else ans=(ans%mod)*qpow(n%mod,mod-2)%mod;
printf("%lld
",ans);
}
return 0;
}