題目：外面有一圈N個結點，中心有一個結點與N個結點都相連，總共就是2*N條邊，刪除N條邊，使N+1個點連通，旋轉相同視為等價，問有多少種情況。
http://acm.hdu.edu.cn/showproblem.php?pid=2481
據說當時現場賽只有清華一個隊過了。非常綜合，其中主要是遞推部分非常難想
好題 ，難！！！！！！！！！！！
做法來源於AC博客：http://hi.baidu.com/aekdycoin/item/517784ec0bf4450b560f1dd1
引用一句話便是：
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位整數的輸入輸出姿勢很是頭疼。。。。。
[cpp]
#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(); www.2cto.com
//構造矩陣
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;
}
作者：ACM_cxlove
