hdu 2865 Birthday Toy 及我对polya的总结

hdu 2865 Birthday Toy 及我对polya的总结

    一直想总结一下这两天学的东西，今天借这个题总结一下。
   正如上篇所说的： 组合计数问题中经常遇到两种困难，第一找出问题通解的表达式，第二是区分讨论问题中哪些应该看成相同的。换句话说，我们就可以将polya 问题分成两部分来分析，从代码上来说，我们也可以分成两部分来分析不同的实现。
    从区分哪些是相同的问题上分析题目，也就是置换群循环节之类的东西。这里分析的时候就不提了。不过就是旋转嘛！
    这个题的重点是如何计数相同珠子不相邻的方案数。
   分析；首先，由于中间一个大圆与每个小圆都相连，所以大圆用去一种颜色之后，只剩下K-1种颜色。
           设K-1种颜色染N个珠子的不同方案数为M，最后就是求M×K mod 1000000007。
           方法跟pku 2888一样，但是这次矩阵的规模很大，所以不能用矩阵来存这个图形了。
           但由于此处规定相邻珠子颜色不同， 则邻接阵为对角线上元素全为0,，其余元素全为1。
           该矩阵的幂的迹，可以推导出公式 ( p - 1 ) ^ n + ( -1 ) ^ n * ( p - 1 )其中p是矩阵的阶数，也就是K-1。
    这个公式是怎么求出来的呢？？？？？？
    几乎所有的日志中都是这个公式，但是没见到有解释怎么求的这个公式，我说说我的想法：
    假设A的n-1次幂为：

其中x_n是对角线上的值。乘以对角线上全0，其余为1的矩阵后。
                                                                                 则1)    x_n = y_n-1*(p-1);
                                                                                    2)    y_n = x_n-1+y_n-1*(p-2);
上面两个式子可以解出来x_n = (p-2)*x_n-1 + (p-1)*x_n-2;
事实上，到这一步就能解了，利用矩阵的乘法，然后快速求出x_n的值，进而求出矩阵幂的迹。
当然到这一步，并没有推导出前面的那个公式，上面的递推公式怎么解呢？
注意到：x_n+x_n-1 = (p-1)*(x_n-1+x_n-2) ;
这样的话就能解出x_n+x_n-1;
接下来解出x_n不是问题了。

至此，我们解决了计数问题的两个难题。

我还想总结一下polya上得几个小的知识点。
一：大数的乘法逆元
        大数的乘法逆元我看到了三种方法

1.  暴力法：  
   

    int i;    for (i=1;;i++) { if (((long long int)(n)*i-an)%M==0) break; }

2. 欧拉函数 利用了欧拉函数
   

   long long inv( long long n ){    return pow( n, M - 2 )%M;}

3. 扩展欧几里得 解同于方程

   //扩展欧几里德void exp_gcd( LL a ,LL b,LL &x,LL &y) {     if( b == 0 ) {         x = 1;         y = 0;     }     else {          exp_gcd( b,a%b,x,y );          LL t;          t = x;          x = y;          y = t - a/b*y;     }}//逆元inline LL getNN(LL x) {        LL now , t;        exp_gcd( x, M,now,t );        return (now%M+M)%M;}

4.今天看到一个求逆元的简洁写法

int64 inv(int64 x) {      //简洁版求逆元      if(x == 1) return 1;      return  inv(MOD%x) * (MOD - MOD/x) % MOD;  }  

原理还没有弄明白，应该是扩展欧几里得吧，我猜

第一种最慢！

这个解题报告已经过去好久了，但是我还是要对这个地方进行补充，众所周知，求a的逆元的时候，a和m互质。但是有过不呢？？如卡特兰数要是取模呢？看一个链接，讲的很好：http://hi.baidu.com/spellbreaker/blog/item/3b04ec27923ed91f8b82a11e.html

对于polya实现时也存在两种方法

1. 循环 时间复杂度O(sqrt(n))

   long long ans = 0;    for (long long i = 1; i*i <= n; i++) if (n % i == 0) { //这个地方可以优化，i*i<=n        ans = (ans + (geter(p,i) * euler(n / i))%M ) % M;        if(i*i != n) ans = (ans + (geter(p,n/i) * euler(i))%M ) % M;    }

   
   
2. 第二种是质因数分解，dfs枚举质因数

   void dfs(int step, int now, int n) {    if (step >= cnt + 1) {        ans = (ans + fun(n % mod, now - 1) * euler(n / now) % mod) % mod;        return;    }    for (int i = 0, t = 1; i <= c[step]; t *= p[step], i++)        dfs(step + 1, now * t, n);}

上面的第二种方法时间复杂度低，但是实现起来也复杂。

最后给出我hdu 2865的代码

#include <iostream>#include <stdio.h>#include <string.h>using namespace std;const long long M = 1000000007;long long n,k;#define PP 10000100long long prime[PP];bool isPrime[PP+1];int size;void getPrime() {    memset(isPrime,true,sizeof(isPrime));    int i;    for(i=2; i<=PP/2; i++) {        if(isPrime[i])  //i是素数            for(int j=i+i; j<=PP; j+=i)                isPrime[j]=0;    }    for(i=2; i<=PP; i++) {        if(isPrime[i])            prime[size++]=i;    }}long long euler(long long n)//求欧拉函数{    long long i,l,t;    l=n;    for (i=0;i<size;i++)    {        t=0;        while (l%prime[i]==0) { t++; l=l/prime[i]; }        if (t>0) n=n/prime[i]*(prime[i]-1);        if (l==1) break;        if (l/prime[i]<prime[i]) { n=n/l*(l-1); break; }    }    return n%M;}long long pow(long long a,long long n){    long long ret=1;    long long A=a;    while(n)    {        if (n & 1)            ret=(ret*A)%M;        A=(A*A)%M;        n>>=1;    }    return ret%M;}long long geter(long long p,long long n){    long long ans = 0;    ans = pow(p-1,n);    if(n%2 == 0) ans = (ans + p-1)%M;    else ans = (ans + M - (p-1) )%M;    return ans;}long long inv( long long n ){    return pow( n, M - 2 )%M;}long long polya(long long p,long long n){    long long ans = 0;    for (long long i = 1; i*i <= n; i++) if (n % i == 0) { //这个地方可以优化，i*i<=n        ans = (ans + (geter(p,i) * euler(n / i))%M ) % M;        if(i*i != n) ans = (ans + (geter(p,n/i) * euler(i))%M ) % M;    }    //乘法逆元！！    return (ans*inv(n))%M;}int main(){    getPrime();    while(scanf("%I64d%I64d",&n,&k) != EOF)        printf("%I64d\n",((long long)k*polya(k-1,n)%M) );    return 0;}