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

最新推荐文章于 2019-01-08 08:27:55 发布

蝶翼的罪

最新推荐文章于 2019-01-08 08:27:55 发布

阅读量505

点赞数

分类专栏：置换

本文链接：https://blog.csdn.net/dieyidezui/article/details/40436913

版权

置换专栏收录该内容

2 篇文章 0 订阅

订阅专栏

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

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

因为有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;
}

//计算polya
LL 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);
    }
}