HDU 2865 Birthday Toy(Polya综合)

题目链接；
HDU 2865 Birthday Toy
题意：
有一个n个珠子的环，中心还有一颗珠子，用k种颜色来染。要求相邻珠子的颜色不同，中心珠子的颜色不能和外围任意一颗珠子的颜色一样，考虑旋转，问本质不同的珠子个数？
数据范围：n≤109,k≤109
分析:
和前面POJ 2888 Magic Bracket类似，这里把限制改成了相邻珠子颜色不同且颜色个数范围变为了109.
对于循环个数k，定义dp[i][0]表示一个循环的的珠子经历i个循环颜色和循环起始珠子颜色一致的方案数，定义dp[i][1]表示颜色不一样(初始化dp[1][0]=1,dp[1][1]=0)，那么可以得到转移方程:

dp[i][0]=dp[i−1][1],

dp[i][1]=dp[i−1][1]∗(m−2)+dp[i][0]∗(m−1)

因为m很大，用矩阵快速幂进行加速。定义矩阵：

A=(0m−11m−2)

B=(10)

C=Ak∗B

那么循环个数为k总的方案数（考虑起始珠子有m种颜色可选）是:ϕ(nk)∗m∗C[1][1]条件是n%k=0,最后还要乘上n模mod的逆元。

#include <iostream>#include <cstdio>#include <cstring>#include <string>#include <algorithm>#include <climits>#include <cmath>#include <ctime>#include <cassert>#define IOS ios_base::sync_with_stdio(0); cin.tie(0);using namespace std;typedef long long ll;const ll mod = (ll)(1e9) + 7;struct Matrix{    int row, col;    ll data[10][10];    Matrix () {}    Matrix (int k) {        row = col = 2;        data[1][1] = 0, data[1][2] = 1;        data[2][1] = k - 2, data[2][2] = k - 3;    }    Matrix operator * (const Matrix& rhs) const {        Matrix res;        res.row = row, res.col = rhs.col;        for(int i = 1; i <= res.row; ++i) {            for (int j = 1; j <= res.col; ++j) {                res.data[i][j] = 0;                for(int k = 1; k <= col; ++k) {                    res.data[i][j] += data[i][k] * rhs.data[k][j];                }                res.data[i][j] %= mod;            }        }        return res;    }    Matrix operator ^ (const int m) const {        Matrix res, tmp;        tmp.row = res.row = row, tmp.col = res.col = col;        memset(res.data, 0, sizeof(res.data));        for(int i = 1; i <= row; ++i) res.data[i][i] = 1;        memcpy(tmp.data, data, sizeof(data));        int mm = m;        while(mm) {            if(mm & 1) res = res * tmp;            tmp = tmp * tmp;            mm >>= 1;        }        return res;    }};inline ll phi(ll x){    ll res = x;    for(ll i = 2; i * i <= x; ++i) {        if(x % i == 0) {            res = res / i * (i - 1);            while(x % i == 0) x /= i;        }    }    if(x > 1) res = res / x * (x - 1);    return res;}ll solve(int len, int k){    Matrix res, tmp;    tmp = Matrix(k);    tmp = tmp ^ len;    res.row = 2, res.col = 1;    res.data[1][1] = 1, res.data[2][1] = 0;    res = tmp * res;    return res.data[1][1];}ll quick_pow(int x, ll m){    ll res = 1, tmp = x;    while(m) {        if(m & 1) res = res * tmp % mod;        tmp = tmp * tmp % mod;        m >>= 1;    }    return res;}int main(){    int n, k;    while(~scanf("%d%d", &n, &k)) {        ll ans = 0;        for(int i = 1; i * i <= n; ++i) {            if(n % i) continue;            ans = (ans + phi(n / i) * solve(i, k) % mod) % mod;            if(n / i == i) continue;            ans = (ans + phi(i) * solve(n / i, k) % mod) % mod;        }        printf("%lld\n",ans * k % mod * (k - 1) % mod * quick_pow(n, mod - 2) % mod);    }    return 0;}