P3307-[SDOI2013]项链【Burnside引理,莫比乌斯反演,特征方程】

1|0正题

题目链接:https://www.luogu.com.cn/problem/P3307


1|1题目大意

n个珠子的一个环形项链,每个珠子有三个1k的整数。

  • 两个珠子不同当且仅当它们不能通过翻转或者旋转得到
  • 两个项链不同当且仅当它们不能通过旋转得到
  • 珠子要求上面的数字互质
  • 项链要求相邻珠子不同

求方案数。


1|2解题思路

珠子的计数和项链计数是没什么关系的,所以两个分开求。

先考虑求珠子有多少种。相当与求有多少种三元组,因为三个数字的排列通过以上操作可以得到任何其他排列。
首先答案包括了三个不同的,有两个相同的,全部都相同的。

也就是设Si表示i元组的数目的话答案就是S3+3S2+2S1(容斥一下可以得到的)

可以先不考虑顺序问题,最后除上一个六就好了。

要求互质所以可以直接上莫反就变为求i=1kμ(i)F(i)。对于每个F(i)来说就是有ki个数字不要求互质的答案。
对于三元组来说F(i)=ki3,二元组就是F(i)=ki2可以和三元组一起解决,一元组就是1
这样整除分块就可以k的处理单次询问了。

现在考虑项链的方案,根据Burnside引理的话就是1ni=1nf(gcd(i,n)),其中f(i)表示大小为i的不动点数量(也就是转i次的话置换中环的数量)。化一下就是1nd|nf(d)φ(nd),需要考虑如何求f

因为是不动点,所以每一个环中的颜色都得相同,然后因为每个环都是相互隔开的,也就是可以视为有gcd(n,i)个珠子头尾相连要涂色,不能有相邻颜色相同。
这个f是有递推式f(i)=f(i1)(m2)+f(i2)(m1)(据说法是第一个是直接在末端插入一个颜色的珠子,第二个是先在末端插入一个和首部相同颜色的珠子再插入一个不同颜色的将它们隔开,这个是第一个无法处理的情况)。

这个可以用矩阵乘法快速递推,但是好像有更简单的方法可以一试,考虑将式子化为f(x)af(x1)=bf(x1)+abf(x2)的形式,这样可以将答案化为一个带幂的函数。

那就有方程{ba=m2ab=m1,然后可以解出a=1或者a=1m。带个1进去做。也就是有f(x)f(x1)=(m1)(f(x1)+f(x2))也就是乘上m1就往后退一步,所以有f(x)f(x1)=(m1)x2(f(1)+f(2))直接带f(1)=0,f(2)=m(m1)进去就有f(x)=m(m1)x1f(x1)
设一个A满足

f(x)+(m1)A=(f(x1)+A)

f(x)+mA=f(x1)

又因为有f(x)=m(m1)x1f(x1)

m(m1)x1+mA=0A=(m1)x1

所以把A带回最初式子,把f(x1)换成f(x)就有式子

f(x)+(m1)x=(f(x)(m1)x1)

f(x)=(1)x(m1)+(m1)x

(这个好像是叫特征方程的玩意)

现在就可以用快速幂算f(x)了,φ比较大,不能直接n求单个,可以先把n质因数分解了然后dfs下去就好了。然后还有因为n很大,有可能是模数P的倍数,这个时候我们求答案模P2,然后除以一个P在和nP搞逆元就好了。

时间复杂度O(k+T(n+m))


1|3code

#include<cstdio> #include<cstring> #include<algorithm> #define ll long long using namespace std; const ll N=1e7+10,XJQ=1e9+7,inv6=833333345000000041ll; ll T,n,m,k,ans,P,cnt,num; ll pri[N/10],mu[N],p[110],c[110]; bool v[N]; void init(){ mu[1]=1; for(ll i=2;i<N;i++){ if(!v[i])pri[++num]=i,mu[i]=-1; for(ll j=1;j<=num&&i*pri[j]<N;j++){ v[i*pri[j]]=1; dans } for(ll i=1;i<N;i++) mu[i]+=mu[i-1]; return; } ll ksc(ll a,ll b,ll p){ a%=p;b%=p; ll tmp=(long double)a*b/p; long double ans=a*b-tmp*p; if(ans>=p)ans-=p; else if(ans<0)ans+=p; return ans; } ll power(ll x,ll b,ll p){ ll ans=1; while(b){ if(b&1)ans=ksc(ans,x,p); x=ksc(x,x,p);b>>=1; } return ans; } void Get_M(ll n){ m=2; for(ll l=1,r;l<=n;l=r+1){ r=n/(n/l);ll k=n/l; (m+=ksc(ksc(ksc(k,k,P),k+3,P),mu[r]-mu[l-1],P))%=P; } m=ksc(m,inv6,P); return; } ll f(ll n) {return power(m-1,n,P)+((n&1)?(P-m+1):(m-1))%P;} void dfs(ll dep,ll x,ll val){ if(dep>cnt){ (ans+=ksc(f(n/x),val,P))%=P; return; } dfs(dep+1,x,val); val*=p[dep]-1;x*=p[dep]; dfs(dep+1,x,val); for(ll i=2;i<=c[dep];i++) val=val*p[dep],x*=p[dep],dfs(dep+1,x,val); return; } signed main() { scanf("%lld",&T); init(); while(T--){ scanf("%lld%lld",&n,&k); if(n%XJQ!=0)P=XJQ; else P=XJQ*XJQ; Get_M(k); // printf("%lld",m); ll x=n;cnt=ans=0; for(ll i=1;i<=num&&pri[i]*pri[i]<=x;i++){ if(x%pri[i]==0){ p[++cnt]=pri[i];c[cnt]=0; while(x%pri[i]==0) c[cnt]++,x/=pri[i]; } } if(x>1)p[++cnt]=x,c[cnt]=1; dfs(1,1,1); if(n%XJQ==0) printf("%lld\n",ans/XJQ*power(n/XJQ,XJQ-2,XJQ)%XJQ); else printf("%lld\n",ans*power(n,XJQ-2,XJQ)%XJQ); } return 0; }

__EOF__

本文作者QuantAsk
本文链接https://www.cnblogs.com/QuantAsk/p/14294579.html
关于博主:退役OIer,GD划水选手
版权声明:本博客所有文章除特别声明外,均采用 BY-NC-SA 许可协议。转载请注明出处!
声援博主:如果您觉得文章对您有帮助,可以点击文章右下角推荐一下。您的鼓励是博主的最大动力!
posted @   QuantAsk  阅读(101)  评论(0编辑  收藏  举报
编辑推荐:
· .NET Core 中如何实现缓存的预热?
· 从 HTTP 原因短语缺失研究 HTTP/2 和 HTTP/3 的设计差异
· AI与.NET技术实操系列:向量存储与相似性搜索在 .NET 中的实现
· 基于Microsoft.Extensions.AI核心库实现RAG应用
· Linux系列:如何用heaptrack跟踪.NET程序的非托管内存泄露
阅读排行:
· TypeScript + Deepseek 打造卜卦网站:技术与玄学的结合
· 阿里巴巴 QwQ-32B真的超越了 DeepSeek R-1吗?
· 【译】Visual Studio 中新的强大生产力特性
· 张高兴的大模型开发实战:(一)使用 Selenium 进行网页爬虫
· 【设计模式】告别冗长if-else语句:使用策略模式优化代码结构
点击右上角即可分享
微信分享提示