HDU 2481 Toy（08成都現場 Polya,遞推，矩陣，數論……）
題目：外面有一圈N個結點，中心有一個結點與N個結點都相連，總共就是2*N條邊，刪除N條邊，使N+1個點連通，旋轉相同視為等價，問有多少種情況。
據說當時現場賽只有清華一個隊過了。非常綜合，其中主要是遞推部分非常難想
好題 ，難！！！！！！！！！！！
引用一句話便是：
PS.此題使用到了:
素數篩選，求解尤拉函式,BurnSide引理,二分模擬乘法,遞推的構造,矩陣二分求冪,置換群,列舉...總之是一個不錯的題目,基本把數論的基本知識全考察了一次.
可想而知。。。。
我們先處理一下有多少種可能，然後再考慮旋轉。將AC部落格裡的遞推整理了一下，重新縷了縷
這裡任意取兩個結點討論a,b。那麼總數便是a,b斷開的種數與a,b連在一起的種數的和。
f(n)表示外圈有n個結點時，而a,b是斷開的種數。
g(n)表示外圈有n個結點時，而a,b是連在一起的種數。
如果a,b之間是斷開的，如果與a直接相連的為k個(加上a自己)，那麼顯然這k個要與其它的保持連通的，與中心必須有一條邊，如果有多條邊就形成環了，顯然不滿足生成樹。另外n-k為f[n-k]種，我們可以列舉k，則f[n]=sigma(i*f[n-i]) (n-1>=i>=0)
如果a,b是連在一起的，如果與a,b相連的為k個(包括a,b)，那麼a,b是相鄰的在這k個位置選擇就有k-1種，而這k個與中心相連的選擇有k種，剩下的與這部分是分開的，則為f[n-k]，所以可以列舉k，最終結果g[n]=sigma(i*(i-1)*f[n-i])
(n-1>=i>=2)
則最終的種數便是T[n]=f[n]+g[n]。
f[n]=sigma(i*f(n-i)) (n-1=>i>=0)
f[n]=f(n-1)+2*f(n-2)+3*f(n-3)……(n-1)*f(1)+n*f[0]
=f(n-1)+f(n-2)+f(n-3)……f(1)+f(0)+(f(n-2)+2*f(n-3)……+(n-1)*f(0))
令s[n]為f[i]的前n項和，則上式可以寫成
f[n]=s[n-1]+f(n-2)+2*f(n-3)……(n-2)f(1)
=s[n-1]+sigma(i*f((n-1)-i)) (n-2=>i>=0)
=s[n-1]+f[n-1] (1) s[n-1]=f[n]-f[n-1]
=s[n-2]+f[n-1]+f[n-1]
=s[n-2]+2*f[n-1]
=f[n-1]-f[n-2]+2*f[n-1] 根據(1)式對s[n-2]變形
=3*f[n-1]-f[n-2] 其中f[0]=1,f[1]=1,f[2]=3,f[3]=8
g(n)=sigma[i(i-1)f(n-i)] (1<=i<n)
=1*2*f[n-2]+2*3*f[n-3]+3*4*f[n-4]……(n-1)*(n-2)*f[1]
則g(n-1)=1*2*f[n-3]+2*3*f[n-4]……+(n-2)*(n-3)*f[1]
則g(n)-g(n-1)=2*f[n-2]+4*f[n-3]……(2*(n-2))*f[1]
=2*(f[n-2]+2*f[n-3]……+f[1])
=2*f[n-1]
這個是最基本的遞推式了。。
g[n]=2*(f[1]+f[2]+f[3]……f[n-1])=2*(s[n-1]-f[0])
=2*(f[n]-f[n-1]-1) 其中f[0]=1
AC引入了f[0]解決了g()的一點小問題，但是他在部落格的推導，寫的時候有一點點問題，如果s[n]包括f[0]那麼g[n]是不等於2*s[n-1]，大神已經完成了重要的推導，這應該是筆誤。
對於f[n]的求法，可以用矩陣快速冪乘解決
{f[n],f[n-1]}={f[1],f[0]}*|3 1|^(n-1)
|-1 0|
而g[n]也就可以順便得到，T[n]就處理完畢了。
然後就是Burnside定理，同樣N比較大，肯定是要用尤拉函式優化，列舉迴圈個數
，這裡不再贅述。
開始的時候覺得MOD在10^9,只要用64位整數，中間部分應該都沒有問題，用了擴充套件歐幾里德求逆元，可是連樣例都過不了，才發現n對MOD是極有可能沒有逆元的，徹底無語了。
只能用(a/b)%c=(a%(b*c))/b。這樣的話把取模就變成了MOD*N，範圍一下子到了10^18，這樣子的話中間的乘法便會溢位64位整數。
所有的大整數相乘都得二分模擬。。。
另外64位整數的輸入輸出姿勢很是頭疼。。。。。
#include<iostream> #include<cstring> #include<queue> #include<cstdio> #include<cmath> #include<algorithm> #include<vector> #include<map> #define N 1000000000 #define inf 1<<29 #define LL long long #define eps 1e-7 #define pb(a) push_back(a) #define ub(v,a) upper_bound(v.begin(),v.end(),a) using namespace std; struct Matrix{ LL m[2][2]; }init; LL MOD; int n; bool flag[40000]={0}; int prime[40000],cnt=0; //由於a,b的範圍都是10^18，二分模擬計算a*b LL MultMod(LL a,LL b){ a%=MOD; b%=MOD; if(b<0) b+=MOD; if(a<0) a+=MOD; LL ret=0; while(b){ if(b&1){ ret+=a; if(ret>=MOD) ret-=MOD; } a=a<<1; if(a>=MOD) a-=MOD; b=b>>1; } return ret; } Matrix operator*(Matrix m1,Matrix m2){ Matrix ans; for(int i=0;i<2;i++) for(int j=0;j<2;j++){ ans.m[i][j]=0; for(int k=0;k<2;k++) ans.m[i][j]=(ans.m[i][j]+MultMod(m1.m[i][k],m2.m[k][j]))%MOD; } return ans; } Matrix operator^(Matrix m1,int b){ Matrix ans; for(int i=0;i<2;i++) for(int j=0;j<2;j++) ans.m[i][j]=(i==j); while(b){ if(b&1) ans=ans*m1; m1=m1*m1; b>>=1; } return ans; } //以上為矩陣快速冪乘 void Prime(){ for(int i=2;i<=sqrt(N+1.0);i++){ if(flag[i]) continue; prime[cnt++]=i; for(int j=2;j*i<=sqrt(N+1.0);j++) flag[i*j]=true; } } int Eular(int n){ int ret=1; for(int i=0;i<cnt&&prime[i]*prime[i]<=n;i++){ if(n%prime[i]==0){ n/=prime[i];ret*=prime[i]-1; while(n%prime[i]==0){n/=prime[i];ret=(ret*prime[i])%MOD;} } } if(n>1) ret*=n-1; return ret%MOD; } //以上為打素數表，求尤拉函式 LL Get_T(int k){ if(k==1) return 1; else if(k==2) return 5; Matrix temp=init^(k-2); LL f=3*temp.m[0][0]+temp.m[1][0]; LL g=2*(f-(3*temp.m[0][1]+temp.m[1][1])-1); return (g+f)%MOD; } //計算T值 LL Polya(){ LL sum=0; int i; //Burnside定理，列舉迴圈個數 for(i=1;i*i<n;i++) if(n%i==0){ sum=(sum+MultMod(Eular(i),Get_T(n/i)))%MOD; sum=(sum+MultMod(Eular(n/i),Get_T(i)))%MOD; } if(i*i==n) sum=(sum+MultMod(Get_T(i),Eular(i)))%MOD; return sum/n; } int main(){ Prime(); //構造矩陣 init.m[0][0]=3;init.m[0][1]=1;init.m[1][0]=-1;init.m[1][1]=0; while(scanf("%d%I64d",&n,&MOD)!=EOF){ MOD=(LL)n*MOD; printf("%I64d\n",Polya()%(MOD/n)); } return 0; }