HDU 5868 ACM/ICPC Regional Dalian Online（polya计数+矩阵快速幂+欧拉函数）

Different Circle Permutation

Time Limit: 3000/1500 MS (Java/Others)    Memory Limit: 262144/262144 K (Java/Others)
Total Submission(s): 367    Accepted Submission(s): 177

Problem Description

You may not know this but it's a fact that Xinghai Square is Asia's largest city square. It is located in Dalian and, of course, a landmark of the city. It's an ideal place for outing any time of the year. And now:
  There are N children from a nearby primary school flying kites with a teacher. When they have a rest at noon, part of them (maybe none) sit around the circle flower beds. The angle between any two of them relative to the center of the circle is always a multiple of2πN but always not 2πN.
  Now, the teacher raises a question: How many different ways there are to arrange students sitting around the flower beds according to the rule stated above. To simplify the problem, every student is seen as the same. And to make the answer looks not so great, the teacher adds another specification: two ways are considered the same if they coincide after rotating.

Input

There are T tests (T≤50). Each test contains one integer N.1≤N≤1000000000 (109). Process till the end of input.

Output

For each test, output the answer mod 1000000007 (109+7) in one line.

Sample Input

4710

Sample Output

3515

Source

2016 ACM/ICPC Asia Regional Dalian Online

        人生第一道polya计数的题目，纪念一下。

        题目大意就是，给你一个n个点的还，然后你可以给每个点涂成黑色或者白色，但是要求任意两个相邻的点不能涂成黑色，问在旋转同构的要求下，本质不同的方案数有多少种。

        看到这种题目，问本质不同的方案数，那么一般就是用polya计数原理了。首先，我们不考虑这个本质不同，我们假设有n个点时，很容易发现方案数f(n)=f(n-1)+f(n-2)。分别表示第n个点颜色为黑色白色，和第n个点为黑色的方案数。那么就成了fibonacci数列了，只不过首项是1和3.。

        接着，我们再考虑旋转异构的形况情况，考虑burnside引理：对于一个置换f，若一个着色方案s经过置换后不变，称s为f的不动点。将f的不动点数目记为C(f)，则可以证明等价类数目为所有C(f)的平均值。考虑这题，总共n个点对应有n个置换，利用burnside引理求出所有C(f)的平均值即可。

        具体来说如何计算这个不动点数目呢？对于第i个置换，我们可以对应有gcd(n,i)个循环，每个循环中相邻两个数字间隔i。例如有6个点第4个置换就有2个循环，0、4、8 mod 6 = 2和1、5、3。然后，我们只需要使每个循环里面所有点颜色相同，那么这样子就可以保证该置换都是不动点。既然我们使每个循环里所有的点颜色相同，那么我们可以把每个循环等价为一个点，如此一来，两个循环的不动点数目就是f(2)。那么总的结果就是sigma(f(gcd(i,n)))的平均数，其中（1<=i<=n）。

        f(n)可以用矩阵快速幂加速，但是递推式还是比较慢。但是看到了gcd，我们就可以继续用数论的知识去优化。

                                                                                                       

        枚举所有的因子d，然后快速幂求f(d)，再求个φ，时间复杂度为O(N^1.5)。具体见代码：

#include<bits/stdc++.h>#define mod 1000000007#define LL long long#define N 3using namespace std;struct matrix{LL a[N][N];void init(){memset(a,0,sizeof(a));}friend matrix operator *(matrix x,matrix y)    {        matrix ans; ans.init();        for(int i=1;i<N;i++)            for(int j=1;j<N;j++)            {                for(int k=1;k<N;k++)                    ans.a[i][j]+=(x.a[i][k]*y.a[k][j])%mod;                ans.a[i][j]%=mod;            }        return ans;    }    friend matrix operator ^(matrix x,LL y)    {        matrix ans;        if (y==0)        {            memset(ans.a,0,sizeof(ans.a));            for(int i=1;i<N;i++) ans.a[i][i]=1;            return ans;        } else while ((y&1)==0) y>>=1,x=x*x;        ans=x; y>>=1;        for(;y!=0;y>>=1)        {            x=x*x; if ((y&1)!=0) ans=ans*x;        }        return ans;    }} x;vector<LL> num;int n;LL fastpow(LL x,LL n)//快速幂{    LL ret=1;    while (n)    {        if (n&1) ret=ret*x%mod;        x=x*x%mod; n>>=1;    }    return ret;}LL phi(LL k)//欧拉函数phi{    LL i,s;    s = k;    for(i = 2;i * i <= k; i++)    {        if(k % i == 0) s = s / i * (i - 1);        while(k % i == 0) k /= i;    }    if(k > 1) s = s / k * (k - 1);    return s;}LL fib(LL n)//fibonacci数列{    if (n==1) return 1LL;    if (n==2) return 3LL;    matrix ans=x^(n-2);    LL res=(ans.a[1][1]*3)%mod+ans.a[1][2];    return res%mod;}int main(){    x.a[1][1]=1;    x.a[1][2]=1;    x.a[2][1]=1;    while(~scanf("%d",&n))    {        num.clear();        if (n==1){puts("2");continue;}        for(int i=1;i*i<=n;i++)        {            if (n%i) continue;            num.push_back(i);            if (i*i!=n) num.push_back(n/i);        }        LL ans=0;        for(int i=0;i<num.size();i++)        {            ans+=fib(num[i])*phi(n/num[i]);            if (ans>=mod) ans%=mod;        }        ans=ans*fastpow(n,mod-2)%mod;//求逆元        printf("%I64d\n",(ans+mod)%mod);    }    return 0;}