HDU 5868 Different Circle Permutation（dp+polya+矩阵快速幂+欧拉函数）

Description
环上均匀分布n个点，每个点可以坐人和不坐人，相邻两点不能同时坐人，旋转后相等视为同一种方案，问有多少种方案
Input
多组用例，每组用例输入一个整数n，以文件尾结束输入
(不超过50组用例，1<=n<=1e9)
Output
对于每组用例，输出合法方案数
Sample Input
4
7
10
Sample Output
3
5
15
Solution
问题为从环上的n个点中选一些点染色，相邻两点不能同时染色
先不考虑旋转同构，用g[i]表示在不考虑旋转同构的前提下给i个点的环染色的合法方案数，先考虑一条链上的情况，用dp[i][0]表示第i个点不染色的方案数，dp[i][1]表示第i个点染色的方案数，那么由以下转移方程：
dp[1][0]=dp[1][1]=1
dp[i][0]=dp[i-1][0]+dp[i-1][1]，dp[i][1]=dp[i-1][0]，i=2,3,…,n
显然可以发现dp[i][0]是斐波那契数列的第i+1项，考虑第n个点是否染色，如果不染色那么可以将这个链连成一个环，方案数dp[n][0]；如果染色，那么第1个点必然不染，而从第2个点开始到第n个点还是可以看做一条链上的染色，方案数dp[n-1][1]=dp[n-2][0]
故有g[n]=f[n-1]+f[n+1]，f[i]表示斐波那契数列第i项，进一步由f[i]=f[i-1]+f[i-2]知g[i]=g[i-1]+g[i-2]，其中g[1]=1,g[2]=3
之后考虑旋转，由Burnside引理知合法方案数

枚举d，矩阵快速幂O(logn)求g，求一个数的欧拉函数也是O(logn)，故总时间复杂度
注意n=1时特判答案是2，因为一个人的时候第1个人和第n个人可以同时染色
Code

#include<cstdio>#include<iostream>#include<cstring>#include<algorithm>using namespace std;typedef long long ll;#define mod 1000000007ll typedef long long ll;struct Mat{    ll mat[3][3];//矩阵     int row,col;//矩阵行列数 };Mat mod_mul(Mat a,Mat b,int p)//矩阵乘法 {    Mat ans;    ans.row=a.row;    ans.col=b.col;    memset(ans.mat,0,sizeof(ans.mat));    for(int i=0;i<ans.row;i++)              for(int k=0;k<a.col;k++)            if(a.mat[i][k])                for(int j=0;j<ans.col;j++)                {                    ans.mat[i][j]+=a.mat[i][k]*b.mat[k][j]%p;                    ans.mat[i][j]%=p;                }    return ans;}Mat mod_pow(Mat a,int k,int p)//矩阵快速幂 {    Mat ans;    ans.row=a.row;    ans.col=a.col;    for(int i=0;i<a.row;i++)        for(int j=0;j<a.col;j++)            ans.mat[i][j]=(i==j);    while(k)    {        if(k&1)ans=mod_mul(ans,a,p);        a=mod_mul(a,a,p);        k>>=1;    }    return ans;}#define maxn 55555int euler[maxn],prime[maxn],res;void Get_euler(int n)//求n以内所有数的欧拉函数值,顺便得到n以内所有素数共res个{    memset(euler,0,sizeof(euler));    euler[1]=1;    res=0;    for(int i=2;i<=n;i++)    {        if(!euler[i])euler[i]=i-1,prime[res++]=i;        for(int j=0;j<res&&prime[j]*i<=n;j++)        {            if(i%prime[j]) euler[prime[j]*i]=euler[i]*(prime[j]-1);            else            {                euler[prime[j]*i]=euler[i]*prime[j];                break;            }        }    }}ll f[maxn];void init(){    Get_euler(50000);    f[1]=1,f[2]=3;    for(int i=3;i<=50000;i++)f[i]=(f[i-1]+f[i-2])%mod;}ll get_f(int n){    if(n<=50000)return f[n];    Mat A;    A.row=A.col=2;    A.mat[0][0]=1,A.mat[0][1]=1,    A.mat[1][0]=1,A.mat[1][1]=0;    A=mod_pow(A,n-2,mod);    return 3ll*A.mat[0][0]+A.mat[0][1];}ll get_euler(int n){    if(n<=50000)return euler[n];    ll ans=n;    for(int i=2;i*i<=n;i++)        if(n%i==0)        {            ans=ans/i*(i-1);            while(n%i==0)n/=i;        }    if(n>1)ans=ans/n*(n-1);    return ans;}ll Mod_pow(ll a,ll b){    a%=mod;    ll ans=1ll;    while(b)    {        if(b&1)ans=ans*a%mod;        a=a*a%mod;        b>>=1;    }    return ans;}int main(){    init();    int n;    while(~scanf("%d",&n))    {        if(n==1)        {            printf("2\n");            continue;        }        ll ans=0;        for(int i=1;i*i<=n;i++)            if(n%i==0)            {                ll temp=(get_f(i))%mod*get_euler(n/i)%mod;                ans=(ans+temp)%mod;                if(i*i!=n)                {                    int j=n/i;                    temp=get_f(j)*get_euler(n/j)%mod;                    ans=(ans+temp)%mod;                }                       }        printf("%I64d\n",ans*Mod_pow(n,mod-2)%mod);    }    return 0;}