HDU 2865 Birthday Toy(Polya综合)

题目链接；
HDU 2865 Birthday Toy
题意：
有一个$n$个珠子的环，中心还有一颗珠子，用$k$种颜色来染。要求相邻珠子的颜色不同，中心珠子的颜色不能和外围任意一颗珠子的颜色一样，考虑旋转，问本质不同的珠子个数？
数据范围：$n\leq 10^9, k \leq 10 ^ 9$
分析:
和前面POJ 2888 Magic Bracket类似，这里把限制改成了相邻珠子颜色不同且颜色个数范围变为了$10^9$.
对于循环个数$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=\begin{pmatrix} 0 & 1 \\ m - 1 & m - 2\end{pmatrix}$$

$$B=\begin{pmatrix} 1 \\ 0\end{pmatrix}$$
$$C=A^{k}*B$$
那么循环个数为 $k$总的方案数（考虑起始珠子有$m$种颜色可选）是: $\phi(\frac{n}{k})*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;
}