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次方就要变成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;
}