HDU-5868 Different Circle Permutation

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 of 2π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
4
7
10

Sample Output
3
5
15

题意为一群孩子按圈排座位。可以看成是给一个项链染色，有两种颜色0和1，同时相邻处不能都染1，若两种染色方案旋转后相同则视为同一种方案。
从题意可以明显的看出，这是一道考察Burnside引理的数论题。
首先，变换方式只有旋转，可以得出，答案应为\( \frac{1}{n}\sum_{i=1}^{n}f(gcd(i,n)) \)。其中f(x)为表示长度为x的一个循环节的染色方案数。

在本题中，考虑最后一个涂1时，则f(x)=f(x-2)，最后一个涂0时，则f(x)=f(x-1)。所以f(x)=f(x-1)+f(x+2)。换言之，斐波那契数列。但是同时必须注意的一点是，f(1)=1。因为在循环节长度为1时，若涂1则项链每个环都是1，不符合题意。所以f(1)=1，f(2)=3，f(3)=4，f(4)=7…本人的做法是将f(0)设为2。同时需要注意的一点是，当n=1是，需要特判。对于斐波那契数列，理所当然的想到要使用矩阵快速幂。

同时，考虑到n的数值非常大，直接使用\( \frac{1}{n}\sum_{i=1}^{n}f(gcd(i,n)) \)计算是不可能的。所以对此式进一步简化为\(\sum _{d \mid n} f(gcd(d,n)) * \phi(n/d)\)，其中\( \phi() \)函数，即欧拉函数直接计算即可。

而对于\( \frac{1}{n} \)，因为10^9+7是质数，直接计算逆元即可。

下面是AC的代码：

  1 #include"cstdio"
  2 #include"cstring"
  3 #include"string"
  4 #include"cstdlib"
  5 #include"vector"
  6 #include"set"
  7 #include"map"
  8 #include"cmath"
  9 using namespace std;
 10 typedef long long LL;
 11
 12 const LL MOD=1000000007;
 13 LL gcd(LL a, LL b)
 14 {
 15     if(b==0) return a;
 16     return gcd(b,a%b);
 17 }
 18 LL powMod(LL x,LL n,LL mod)
 19 {
 20     LL res=1;
 21     while(n>0)
 22     {
 23         if(n&1) res=res*x % mod;
 24         x=x*x % mod;
 25         n>>=1;
 26     }
 27     return res;
 28 }
 29 const LL MAXR=100;
 30 const LL MAXC=100;
 31 struct Matrix
 32 {
 33     LL m[MAXR][MAXC];
 34     LL r,c;
 35     Matrix()
 36     {
 37         r=0,c=0;
 38         memset(m,0,sizeof(m));
 39     }
 40 };
 41 //若不能相乘，返回空矩阵
 42 Matrix operator * (const Matrix &a,const Matrix &b)
 43 {
 44     Matrix ans;
 45     LL r=a.r,c=b.c,x=0;
 46     if(a.c!=b.r) return Matrix();
 47     x=a.c;
 48     ans.r=a.r;
 49     ans.c=b.c;
 50     for (LL i=1; i<=r; i++)
 51         for (LL j=1; j<=c; j++)
 52         {
 53             ans.m[i][j]=0;
 54             for (LL k=1; k<=x; k++)
 55                 ans.m[i][j]=(ans.m[i][j]+a.m[i][k]*b.m[k][j]) % MOD;
 56         }
 57     return ans;
 58 }
 59 //若不能乘方，返回空矩阵
 60 Matrix operator ^ (Matrix &base,LL pow)
 61 {
 62     if(base.r!=base.c) return Matrix();
 63     LL x=base.r;
 64     Matrix ans;
 65     ans.r=ans.c=x;
 66     for (LL i=1; i<=x; i++) ans.m[i][i]=1;
 67     while (pow)
 68     {
 69         if (pow&1) ans=ans*base;
 70         base=base*base;
 71         pow>>=1;
 72     }
 73     return ans;
 74 }
 75 LL eularPhi(LL n)
 76 {
 77     LL res=n;
 78     for(LL i=2;i*i<=n;i++)
 79     {
 80         if(n%i==0)
 81         {
 82             res=res/i*(i-1);
 83             for(;n%i==0;n/=i);
 84         }
 85     }
 86     if(n!=1) res=res/n*(n-1);
 87     return res;
 88 }
 89 //拓展欧几里得算法
 90 //求ax+by=gcd(a,b)的一组解
 91 //其他解为x=x0+kb‘,y=y0-ka‘
 92 //a‘=a/gcd(a,b),b‘=b/gcd(a,b)
 93 LL extgcd(LL a,LL b,LL &x,LL &y)
 94 {
 95     LL d=a;
 96     if(b!=0)
 97     {
 98         d=extgcd(b,a%b,y,x);
 99         y-=(a/b)*x;
100     }
101     else { x=1; y=0; }
102     return d;
103 }
104 //求逆元
105 //a和ｍ应该互质
106 LL modInverse(LL a,LL m)
107 {
108     LL x,y;
109     extgcd(a,m,x,y);
110     return (m+x%m)%m;
111 }
112 LL f(LL x)
113 {
114     Matrix a,b;
115     a.r=1; a.c=2; a.m[1][1]=2; a.m[1][2]=1;
116     b.r=2; b.c=2; b.m[1][1]=0; b.m[1][2]=b.m[2][1]=b.m[2][2]=1;
117     return (a*(b^x)).m[1][1];
118 }
119 int main()
120 {
121 #ifdef LOCAL
122     freopen("in.txt","r",stdin);
123 #endif
124     LL n;
125     while(scanf("%lld",&n)!=EOF)
126     {
127         if(n==1)
128         {
129             printf("%d\n",2);
130             continue;
131         }
132         LL ans=0;
133         for(LL i=1;i<=sqrt(n);i++)
134             if(n%i==0)
135             {
136                 LL d=i;
137                 ans=(ans+f(d)*eularPhi(n/d))%MOD;
138                 d=n/i;
139                 if(d*d==n) continue;
140                 ans=(ans+f(d)*eularPhi(n/d))%MOD;
141             }
142         printf("%lld\n",ans*modInverse(n,MOD)%MOD);
143     }
144     return 0;
145 }