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