【BZOJ3202】項鏈（莫比烏斯反演，Burnside引理）
【BZOJ3202】項鏈（莫比烏斯反演，Burnside引理）
題面
BZOJ
洛谷
題解
首先讀完題目，很明顯的感覺就是，分成了兩個部分計算。
首先計算本質不同的珠子個數，再計算本質不同的項鏈個數。
前面一個部分和\(gcd\)相關，一種莫比烏斯反演的感覺。
後面一個部分出現了旋轉操作，要求本質不同的方案數，不難想到Burnside引理。
首先先考慮怎麽求本質不同的珠子個數。
我們直接考慮無序的三元組\((x,y,z)\)，滿足\(x,y,z\le a,gcd(x,y,z)=1\)
容斥考慮最終答案，就是\(\frac{1}{6}(S3-S2*3+S1*2)\)，其中\(S3\)是無序的三元組個數，\(S2\)
\(S1\)顯然只有一個，即可以構成唯一合法三元組\((1,1,1)\)
考慮怎麽計算\(S2\)和\(S3\)。
看起來\(s2\)很熟悉，也就是\(\sum_{i=1}^a\sum_{j=1}^a[gcd(i,j)=1]\)
然後直接莫比烏斯反演。
設\(F(x)=\sum_{i=1}^a\sum_{j=1}^a[gcd(i,j)=x]\)
\(G(x)=\sum_{x|d}F(d)=\sum_{i=1}^a\sum_{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}]^2\mu(i)\)。
然後是\(S3\)，即三元組個數，和上面一樣的推法，可以得到貢獻是\(\sum_{i=1}^a[\frac{a}{i}]^3\mu(i)\)
好了，那麽計算本質不同的珠子的個數就算完了。
接下來考慮如何計算本質不同的項鏈的個數。
因為只有旋轉計算本質不同，那麽置換一共\(n\)個，是順時針旋轉\(i\)位。
那麽答案就是所有置換的不動點個數除以\(n\)。那麽不難往Polya定理方面靠。
只需要知道每個循環的大小即可。
對於順時針旋轉\(i\)
因為是不動點的個數，因此每一個循環內的所有珠子的顏色必須相同。同時，題目限制相鄰兩個位置的珠子顏色不能相同。因此，這裏可以等價的看做有\(gcd(n,i)\)個珠子，收尾相連，要求相鄰的珠子不同色的方案數，假裝這個的方案數是\(f(x)\),表示\(x\)個珠子收尾相連時候的方案數。那麽現在要求的答案就是\(\sum_{i=1}^nf(gcd(i,n))\)。\(n\)太大了，考慮優化。
我們枚舉\(gcd\),即\(\sum_{d=1}^n\sum_{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)\)
那麽列出方程組：
\[\begin{cases}b-a&=m-2\\ab&=m-1\end{cases}\]
所以有\(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)\)
一般這樣子都可以構建一個遞推式。
\[
\begin{aligned}
f(x)&=m(m-1)^{x-1}-f(x-1)\f(x)+A(m-1)&=-(f(x-1)+A)\mA&=m(m-1)^{x-1}\A&=(m-1)^{x-1}
\end{aligned}
\]
所以，我們可以構建出這樣一個等式：
\[f(x)-(m-1)^{x}=-(f(x)-(m-1)^{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\n",ans);
}
return 0;
}【BZOJ3202】項鏈（莫比烏斯反演，Burnside引理）