HDU 2481 Toy（08成都现场 Polya,递推，矩阵，数论……）
题目：外面有一圈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;}
