枚舉一個長度爲\(n\)爲迴文串,它的全部循環位移均可以產生貢獻。
可是這樣算重了。重複的地方在於可能多個迴文串循環同構,或者可能有的迴文串通過小於\(n\)次循環位移後可以獲得自身。
一個比較好的處理方式是:對每一個迴文串求最小的\(x\)使這個串通過\(x\)次循環位移後能夠再次成爲一個迴文串。這樣對每一個迴文串求\(\sum x\)顯然就不會算重了。
考慮一個串的\(x\)是什麼。顯然會和這個串的最小循環節長度有關。實際上若是最小循環節長度爲偶數,那麼\(x\)就會是這個長度的一半;不然就等於這個長度。算法
形式化地,若是一個迴文串的最小循環節長度爲\(i\),那麼它對答案的貢獻就是\(h(i)=i\frac{1+[i\mbox{是奇數}]}{2}\)。函數
設最小循環節爲\(i\)的迴文串共有\(f(i)\)個,那麼咱們要求的答案就是spa
又由於$$\sum_{d|n}f(d)=k^{\lceil\frac n2\rceil}=g(n)$$code
因此$$f(n)=\sum_{d|n}g(d)\mu(\frac nd)$$get
代入原式$$Ans=\sum_{d|n}\sum_{i|d}g(i)\mu(\frac di)h(d)\=\sum_{i|n}g(i)\sum_{d|\frac ni}\mu(d)h(id)$$it
咱們但願能夠把\(h(id)\)中的\(i\)提出來,這樣後半部分就是一個關於\(\frac ni\)的函數了。io
由於\(h(x)\)不是\(x\)就是\(\frac x2\),咱們發現\(h(id)\neq d\times h(i)\)當且僅當\(i\)是奇數且\(d\)是偶數,而\(d|\frac ni\)因此\(d\)是偶數就說明\(\frac ni\)也是偶數。那麼咱們如今假設\(i\)是奇數且\(\frac ni\)是偶數,考慮下面這個式子的取值。ast
顯然只有\(\mu(d)\)非零項有貢獻,而\(\frac ni\)中含有\(2\)這個因子就使得全部\(\mu(d)\)非零項中含\(2\)與不含\(2\)的\(d\)能夠一一對應。他們的\(h(id)\)的值是相同的,而\(\mu(d)\)的值剛好相反,因此這個式子的值必定爲\(0\)。class
話說回來。咱們如今已經知道了\(h(id)\neq d\times h(i)\)的狀況沒有貢獻,就能夠放心大膽地用這一種變換了。
(在枚舉\(i\)時須要跳過\(i\)是奇數而\(\frac ni\)是偶數的項)
考慮後面的東西是個啥。仍是隻有\(\mu(d)\)非零項有貢獻,也就是說含有奇數個質因子的\(d\)會乘上\(-1\),含有偶數個質因子的\(d\)爲乘上\(1\),因此這個值至關因而將\(\frac ni\)質因數分解爲\(p_1^{a_1}p_2^{a_2}...p_k^{a_k}\)後,爲\(\prod_{i=1}^k(1-p_i)\)。
因此到這裏就比較簡單了。先用\(Pollard-Rho\)算法將\(n\)分解,再dfs枚舉\(n\)的每個約數\(d\),在搜索的過程當中天然能夠求出那個\(\prod_{i=1}^k(1-p_i)\)。
#include<cstdio> #include<algorithm> #include<ctime> using namespace std; #define ll long long ll mul(ll x,ll y,ll m){ x%=m;y%=m; return (x*y-(ll)(((long double)x*y+0.5)/(long double)m)*m+m)%m; } ll fastpow(ll x,ll y,ll m){ ll res=1; while (y) {if (y&1) res=mul(res,x,m);x=mul(x,x,m);y>>=1;} return res; } ll f[]={2,3,5,7,11,13,17,19,23,29}; bool MR(ll p){ for (int i=0;i<10;++i){ if (p<=f[i]) break; if (fastpow(f[i],p-1,p)!=1) return false; ll pp=p-1; while (~pp&1){ pp>>=1;ll y=fastpow(f[i],pp,p); if (mul(y,y,p)==1&&y!=1&&y!=p-1) return false; } } return true; } ll PR(ll n,ll c){ ll i=0,k=2,x,y;x=y=1+rand()%(n-1); while (1){ x=(mul(x,x,n)+c)%n; ll d=__gcd((y-x+n)%n,n); if (d!=1&&d!=n) return d; if (x==y) return n; if (++i==k) y=x,k<<=1; } } ll tmp[100];int len; void fact(ll n){ if (n==1) return; if (MR(n)) {tmp[++len]=n;return;} ll p=n;for (int c=233;p==n;--c) p=PR(p,c); fact(p);fact(n/p); } int Case,q[100],cnt,mod;ll n,k,p[100],ans; int fpow(int x,ll y){ ll res=1; while (y) {if (y&1) res=1ll*res*x%mod;x=1ll*x*x%mod;y>>=1;} return res; } int g(ll n){return fpow(k,(n+1)>>1);} int h(ll n){return (n&1?n:n>>1)%mod;} void dfs(int i,ll d,int pro){ if (i==cnt+1){ if ((n/d&1)&&(d&1)==0) return; ans=(ans+1ll*g(n/d)*h(n/d)%mod*pro)%mod; return; } dfs(i+1,d,pro);pro=1ll*pro*(mod+1-p[i]%mod)%mod; for (int j=1;j<=q[i];++j) d*=p[i],dfs(i+1,d,pro); } int main(){ srand(20020415); scanf("%d",&Case);while (Case--){ scanf("%lld%lld%d",&n,&k,&mod);k%=mod; len=cnt=0;fact(n);sort(tmp+1,tmp+len+1); for (int i=1;i<=len;++i){ if (tmp[i]!=tmp[i-1]) p[++cnt]=tmp[i],q[cnt]=0; ++q[cnt]; } ans=0;dfs(1,1,1);printf("%lld\n",ans); } return 0; }