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;
}