HDU 2481 Toy(Burnide引数)

テーマリンク:http://acm.hdu.edu.cn/showproblem.php?pid=2481
外にはNの結点があります.中心にはNの結点と繋がっています.全部で2*Nの辺です.Nの辺を削除して、N＋1の点を通して、回転は同じです.等価として、どれぐらいの状況がありますか?
考え方:ソースアドレス:http://blog.csdn.net/acm_cxlove/articale/detail/7868589 まず外輪のn個の点がどれぐらいの可能性があるかを計算してから回転を考えます.ここで任意に二つの結点a、bを取って討論します.では、総数はaであり、bは切断された種数とa、bは結合された種数の和である(この分析方法は重要である).f(n)は外輪にn個の結点がある場合を示し、a,bは断線の種数である.g(n)は外輪にn個の結点がある場合を示し、a,bは連係する種数である.
(1)a,bの間が切断されている場合(下の左側の図のように)、aと直接に接続されているのがk個(a自身を加えて)であると、明らかにこのk個は他の保持と連結されているものであり、中心とは一つの辺が必要であり、複数の辺があればリングが形成されており、明らかに生成木を満足しない.なお、n-kはf[n-k]種であり、このn-k種は断続的であることは間違いない.kを列挙すると、f[n]=sigma(k＊f[n-k])(1<=k==n)となります.

(2)a,bが接続されている場合(上の右側の図のように)、a,bに接続されているものがk個(a,bを含む)であれば、a,bは隣接しているものがk-1種あり、このk個は中心に接続されているものがk種あり、残りはこの部分と分離されているものがf[n-k]であるので、kを列挙しても良い.
最終的な数はT[n]=f[n]+g[n]である.
f[n]=sigma(k＊f[n-k])(1<=k<=n)から:
f[n]=1*f(n-1)+2*f(n-2)+3*f(n-3)+…+(n-1)*f(1)+n*f[0]
    =s[n-1]+1*f(n-2)+2*f(n-3)++(n-2)*f(1)+(n-1)*f[0]
    =s[n-1]+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]     =3＊f[n-1]-f[n-2]   ここで、f[0]=1,f[1]=1,f[2]=3,f[3]=8 g[n] =sigama((k-1)*k*f[n-k](2<=k<=n) ここからg[1]=0が得られます.      =1*2*f[n-2]+2*3*f[n-3]+3*4*f[n-4]+…+(n-1)*(n-2)*f[1]+n*(n-1)*f[0]g[n-1]=          1*2*f[n-3]+2*3*f[n-4]+…+(n-2)*(n-3)*f[1]+(n-1)*(n-1)*f[0]だから:g(n)-g(n-1)=2*f[n-2]+4*f[n-3]…(2*(n-2)*f[1]*f[1]……(n-1)*                 =2*(1＊f[n-2]+2*f[n-3]+…+(n-2)*f[1]+(n-1)*f[0]                 =2＊f[n-1]をさらに入手すると、g[n]=2*(f[1]+f[2]+f[3]......f[n-1]+g[1]=2*(s[n-1]-f[0]+g[1]が得られます.              =2*(f[n]-f[n-1]-1)+g[1]   
              =2*(f[n]-f[n-1]-1) f[0]=1、g[1]=0はf[n]の求法に対して、行列で素早くべき乗して解決できます.マトリクスAを定義します.

f[n],f[n-1]=(f[1],f[0]*A^(n-1)となります.g[n]もついでに入手できます.T[n]は処理済みです.
そしてBurnideの定理です.同じNの方が大きいです.必ずEurer関数を使って最適化し、エニュメレート・サイクルの数を指定します.n対MODは逆元がない可能性があります.(a/b)%c=(a%(b*)/bを使います.このように型取りをMOD*Nに変えて10^18までの範囲にすると、中間の掛け算は64桁の整数が溢れます.二分をかける

const int N=2;

const int M=40000;

int n;

i64 mod;



class Matrix

{

public:

    i64 a[N][N];



    Matrix()

    {

        clr(a,0);

    }



    void init(int x)

    {

        if(x==0) a[0][0]=a[0][1]=a[1][0]=a[1][1]=0;

        else if(x==1) a[0][0]=a[1][1]=1,a[0][1]=a[1][0]=0;

    }

};



i64 mul(i64 a,i64 b)

{

    a=(a%mod+mod)%mod;

    b=(b%mod+mod)%mod;

    i64 ans=0;

    while(b)

    {

        if(b&1) ans=(ans+a)%mod;

        a=(a+a)%mod;

        b>>=1;

    }

    return ans;

}



Matrix operator*(Matrix &a,Matrix &b)

{

    int i,j,k;

    Matrix p;

    p.init(0);

    FOR0(i,2) FOR0(j,2) FOR0(k,2)

    {

        (p.a[i][j]+=mul(a.a[i][k],b.a[k][j])%mod)%=mod;

    }

    return p;

}



Matrix operator^(Matrix a,int t)

{

    Matrix ans;

    ans.init(1);

    while(t)

    {

        if(t&1) ans=ans*a;

        a=a*a;

        t>>=1;

    }

    return ans;

}



Matrix a,b;

int prime[M],tag[M],cnt;

vector<int> factor;





void init()

{

    int i,j;

    for(i=2;i<180;i++) if(!tag[i])

    {

        for(j=i*i;j<M;j+=i) tag[j]=1;

    }

    for(i=2;i<M;i++) if(!tag[i])

    {

        prime[++cnt]=i;

    }

}



int eular(int n)

{

    int ans=n,i;

    for(i=1;i<=cnt&&(i64)prime[i]*prime[i]<=n;i++) if(n%prime[i]==0)

    {

        ans=ans/prime[i]*(prime[i]-1);

        while(n%prime[i]==0) n/=prime[i];

    }

    if(n>1) ans=ans/n*(n-1);

    return ans;

}



i64 get(int n)

{

    if(n==1) return 1;

    if(n==2) return 5;

    a=b^(n-2);

    i64 f=3*a.a[0][0]+a.a[1][0];

    i64 g=2*(f-(3*a.a[0][1]+a.a[1][1])-1);

    return (g+f)%mod;

}



int main()

{

    init();

    b.a[0][0]=3; b.a[0][1]=1; b.a[1][0]=-1; b.a[1][1]=0;

    while(scanf("%d%I64d",&n,&mod)!=-1)

    {

        mod*=n;

        i64 ans=0;

        int i;

        for(i=1;i*i<=n;i++) if(n%i==0)

        {

            (ans+=mul(eular(i),get(n/i))%mod)%=mod;

            if(i*i==n) continue;

            (ans+=mul(eular(n/i),get(i))%mod)%=mod;

        }

        ans=ans/n%(mod/n);

        printf("%I64d
",ans);

    }

    return 0;

}