hdu 5868 Different Circle Permutation Polya定理　欧拉函数优化

题意：叉姐给出的解释很明白，就是n个点围成一个圈，每个人都可染成黑或者白，任意相邻的两个人不可以染成黑色，并且循环同构，问染色的方案数

思路：如果这道题实在弄不懂，可以先按顺序把下面两题AC了

http://blog.csdn.net/wyt734933289/article/details/52472495

http://blog.csdn.net/wyt734933289/article/details/52540674

把上面两道题完成后，基本上思路方向是一致了，但还要其他分析

我们考虑下普通polya定理是怎么做的，先找出一种不考虑旋转同构的方案数，找到后把方案数乘上此等价类里面有几个元素，这题也一样，设f(n)为不考虑旋转同构时的方案数，f(1) = 1, f(2) = 3, f(3) = 4, f(4) = 7(这里有个题外话，f(1)应该是等于1的，因为此时2pi/1，但是输出的时候又要特判n == 1输出2，这个小细节知道就好了，做题的思路最重要)，可以看到f(n) = f(n - 1) + f(n - 2)  (n >= 3)，可以看到和斐波那契数列形式一样，那么就可以用矩阵快速幂加速，这里有个地方要注意下斐波那契数列的形式是这样的

(抱歉，不太会画这个矩阵)

fn+2        (1    0 )         (f1)

           =                 ^n

fn+1        (1    0)           (f0)

而下面那个举着是从f[1]开始的，就是f[0]的位置有f[1]代，f[1]由f2[]代，自然前面就会变成要想f[n+2]变成f
，那么n次方就要变成n - 2次方

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=5868

#include <cstdio>
#include <cstring>
#include <cmath>
#include <vector>
#include <algorithm>
using namespace std;

typedef long long ll;
const ll mod = 1000000007;

struct Matrix
{
ll mat[5][5];
int r, c;

void init(int r, int c)
{
this->r = r, this->c = c;
memset(mat, 0, sizeof(mat));
}

void identity()//化成单位矩阵
{
for(int i = 0; i < r; i++)
mat[i][i] = 1;
}
};

Matrix operator * (const Matrix &a, const Matrix &b)
{
Matrix c;
c.init(a.r, b.c);
for(int i = 0; i < a.r; i++)
for(int k = 0; k < a.c; k++)
if(a.mat[i][k])//矩阵乘法优化
for(int j = 0; j < b.c; j++)
c.mat[i][j] = (c.mat[i][j] + a.mat[i][k] * b.mat[k][j] % mod) % mod;
return c;
}

Matrix operator ^ (Matrix x, int n)
{
Matrix res;
res.init(x.r, x.c), res.identity();
while(n)
{
if(n & 1) res = res * x;
x = x * x;
n >>= 1;
}
return res;
}

vector<ll> divisors(ll n)
{
vector<ll> res;
ll m = (ll)floor(sqrt(n * 1.0) + 0.5);
for(ll i = 1; i <= m; i++)
{
if(n % i == 0)
{
res.push_back(i);
if(i != n / i) res.push_back(n / i);
}
}
return res;
}

vector<ll> prime_factor(ll n)
{
vector<ll> res;
ll m = (ll)floor(sqrt(n * 1.0) + 0.5);
for(ll i = 2; i <= m; i++)
{
if(n % i == 0)
{
res.push_back(i);
while(n % i == 0) n /= i;
}
}
if(n > 1) res.push_back(n);
return res;
}

ll euler_phi(ll n, const vector<ll> &prime)
{
ll res = n;
int m = prime.size();
for(int i = 0; i < m; i++)
{
if(n % prime[i] == 0)
{
res = res / prime[i] * (prime[i] - 1);
while(n % prime[i] == 0) n /= prime[i];
}
}
return res % mod;
}

ll quick_pow(ll x, ll n)
{
ll res = 1;
x %= mod;
while(n)
{
if(n & 1) res = res * x % mod;
x = x * x % mod;
n >>= 1;
}
return res;
}

ll f(ll n)
{
if(n == 1) return 1;
else if(n == 2) return 3;
else if(n == 3) return 4;
else if(n == 4) return 7;
else
{
Matrix a, b, c;
a.init(2, 2), b.init(2, 1), c.init(2, 1);
a.mat[0][0] = a.mat[0][1] = a.mat[1][0] = 1;
b.mat[0][0] = 3, b.mat[1][0] = 1;
n -= 2;//这里可以－１，那答案就是c.mat[1][0] % mod
a = a ^ n;
c = a * b;
return c.mat[0][0] % mod;
}
}

ll polya(ll n)
{
vector<ll> divs = divisors(n);
vector<ll> prime = prime_factor(n);
int m = divs.size();
ll ans = 0;
for(int i = 0; i < m; i++)
{
ll euler = euler_phi(divs[i], prime);
ans += euler * f(n / divs[i]) % mod;
ans %= mod;
}
return ans * quick_pow(n, mod - 2) % mod;//根据费马小定理，这里相当于求n在mod1000000007下的逆元
}

int main()
{
ll n;
while(~scanf("%lld", &n))
{
if(n == 1)
{
printf("2\n");
continue;
}
ll ans = polya(n);
printf("%lld\n", ans);
}
return 0;
}