原题链接为：HDU-5868
题意为一群孩子按圈排座位。可以看成是给一个项链染色，有两种颜色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
",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
",ans*modInverse(n,MOD)%MOD);
143 }
144 return 0;
145 }