HDU 2865 Birthday Toy polya 矩阵快速幂 欧拉函数

分析：题意略过。就是求相邻颜色不同的置换数目。还是老思路，求不动点然后取平均数。下面讲如何求不动点。

假设我们知道了一个置换的循环数目。如何求相邻颜色不同的不动点数目。我们把环从正上方切开，然后取正上方的那个珠子加到尾部。这个问题的本质是一样的。

因为有k种颜色，内圈的珠子取一个，外圈只能取剩下的k-1中，下面我们就认为k = k-1,不考虑内圈。

 从左往右推，我们假设第i个珠子的取第j种颜色的方案数为c(i,j),那么有这样的公式  c(i,j) =  ∑c(i-1,p)( (p = 1,2,...k,  p != j);根据这个我们可以推出一个矩阵

0 1 1       1

1 0 1       0

1 1 0       0

右边那列是第一个珠子的取值情况，然后每次左乘一个矩阵，就变成下一个珠子的取值情况，而方案数，就是最后一个珠子的取值的总和。但是这样分析后我们发现这样是不够的，因为矩阵的幂  =  循环数 最大为n， 矩阵的长 = K 颜色数，而这两者取值都可以到1e9。如果这样解的话，时间复杂度为O(n^3*log k);我们给行列一些标号然后再研究他们的规律且看

设矩阵长为 n

0 1 1     xi      xi+1 = yi * (n-1)

1 0 1  * yi  =  yi+1 = xi + yi *(n-2)

1 1 0     yi      yi+1 = xi + yi *(n-2)

 所以，与其计算之前那个矩阵的n-1次幂然后再乘以右边那一列，不如直接计算xi yi;回到之前， 当我们计算矩阵的n-1次幂然后乘以那一列后，除了第一行的数，其他的加起来就是结果。换言之，y(n-1)*(n-1) 正好等于xn，我们计算新的矩阵 

0 n-1

1 n -2 

的n次幂直接取第一行第一列的数即可。

解决了这个问题，我们继续看，如果用常规枚举GCD的方法，需要O（n）的时间复杂度，也是无法承受，所以我们要利用欧拉函数来求1-n的所有GCD，详见代码。

#include<cstdio>typedef long long LL;const int SIZE = 2;const int MOD = 1000000007;const int N = 40009;const int M = 1000008;LL prime[N];int phi[M];bool vi[N];int p;//数论相关void get_prim_and_eular(){    p = 0;    for(int i = 2; i < N; i++){        if(!vi[i]){            prime[p++] = i;            for(int j = i * i; j < N; j += i) vi[j] = true;        }    }    phi[1] = 1;    for(int i = 2; i < M; i++) if(!phi[i]){        for(int j = i; j < M; j += i){            if(!phi[j]) phi[j] = j;            phi[j] = phi[j]/ i * (i-1);        }    }}int eular(int n){    if(n < M) return phi[n];    int ans = n;    for(int i = 0; i < p && prime[i] * prime[i] <= n; i++) if(n % prime[i] == 0){        ans = ans / prime[i] * (prime[i]-1);        while(n % prime[i] == 0) n /= prime[i];    }    if(n > 1) ans = ans / n * (n-1);    return ans;}LL pow_mod(LL base, int exp){    LL ret = 1;    while(exp){        if(exp&1) ret = ret * base % MOD;        base = base * base % MOD;        exp >>= 1;    }    return ret;}LL inver(int t){ return pow_mod(t, MOD-2);}//快速幂相关struct Mat{    LL m[SIZE][SIZE];    Mat(){m[0][0] = m[0][1] = m[1][0] = m[1][1] = 0;}};Mat mul(const Mat &A,const Mat &B){    Mat C;    for(int i = 0; i < SIZE; i++){        for(int j = 0; j < SIZE;j++){            for(int k = 0; k < SIZE; k++)                C.m[i][j] = (C.m[i][j] + A.m[i][k] * B.m[k][j]) % MOD;        }    }    return C;}Mat pow(Mat A, LL n){    Mat B;    for(int i = 0; i < SIZE; i++) B.m[i][i] = 1;    while(n){       if(n&1) B = mul(B, A);       A = mul(A, A);       n >>= 1;    }    return B;}//计算polyaLL cal(int d, int c){    Mat a;    a.m[0][1] = c - 1;    a.m[1][0] = 1;    a.m[1][1] = c - 2;    return c * pow(a, d).m[0][0] % MOD;}int main(){    get_prim_and_eular();    int n,k;    while(scanf("%d%d", &n, &k) == 2){        LL ans = 0;        int t = n;        for(int i = 1; i * i<= n; i++)if(n % i == 0){            ans = (ans + cal(i, k-1) * eular(n/i) % MOD) % MOD;            if(i * i != n) ans = (ans + cal(n/i, k-1) * eular(i) % MOD) % MOD;        }        ans = ans * inver(n) % MOD;        ans = ans * k % MOD;        printf("%I64d\n", ans);    }}