HDU Birthday Toy 特殊限制的polya

題意：用k中顏色給n個珠子塗色，注意相鄰的兩個珠子不能途同樣的顏色，求不等價著色數mod 10000000007。另外這道題還有另一個條件，n個小珠子一定圍著一個大珠子，所以必須從k中顏色裡面選擇一種顏色給大珠子著色。

由於此題資料較大1<= k, n <= 10^9，顯然不能像POJ2888那樣用矩陣求迴路，但是為了便於分析，我們還是先從矩陣入手。由於題目本身的特點：相鄰的兩個珠子不能同色。

那麼主對角線上的元素全為0，其餘全為1。而我們要的就是矩陣求冪後主對角元素的值。····分析很簡單，但是畫出來比較麻煩，直接借用：http://blog.csdn.net/wukonwukon/article/details/7215467童鞋的：

分析；首先，由於中間一個大圓與每個小圓都相連，所以大圓用去一種顏色之後，只剩下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) ;  （等式1）
這樣的話就能解出x_n+x_n-1;
接下來解出x_n不是問題了。
補充：由等式1可以得到(x_2+x_1) = p-1 (這是顯然的，應為x_1 = 0)
(x_3+x_2) = (p-1)^2
(x_4+x_3) = (p-1)^3
····(x_n+x_n-1) = (p-1)^(n-1)迭代一下x_n = (p-1)^(n-1) - (p-1)^(n-2) + (p-1)^(n-3) - (p-1)^(n-4) ·····-(p-1)^2
+ (p-1)OK，等比數列求和。所以( p - 1 ) ^ n + ( -1 ) ^ n * ( p - 1 ) 是所有主對角元素之和。

#include<cmath>#include<cstring>#include<algorithm>#include<iostream>#include<cstdio>using namespace std;#define lint __int64const int MAXN = 1000009;const int modulo = 1000000007;int a[MAXN], p[MAXN], pn;int eul[MAXN];lint n, k, ans;void Prime(){    int i, j; pn = 0;    memset(a,0,sizeof(a));    for ( i = 2; i < MAXN; i++ )    {        if ( !a[i] ) p[pn++] = i;        for(j = 0; j < pn && i*p[j] < MAXN && (p[j]<=a[i] || a[i]==0); j++)            a[i*p[j]] = p[j];    }}void Euler(){    eul[1] = 1;    for ( int i = 2; i < MAXN; i++ )    {        if ( a[i] == 0 )            eul[i] = i - 1;        else        {            lint k = i / a[i];            if ( k % a[i] == 0 ) eul[i] = eul[k] * a[i];            else eul[i] = eul[k] * (a[i]-1);        }    }}int Euler ( int n ){    if ( n < MAXN )        return eul[n] % modulo;    int i, ret = n;    for ( i = 0; i < pn && p[i] * p[i] <= n; i++ )    {        if ( n % p[i] == 0 )        {            ret = ret - ret / p[i];            while ( n % p[i] == 0 ) n /= p[i];        }    }    if ( n > 1 )        ret = ret - ret / n;    return ret % modulo;}lint mod_exp ( lint a, lint b ){    lint ret = 1;    a = a % modulo;    while ( b >= 1 )    {        if ( b & 1 )            ret = ret * a % modulo;        a = a * a % modulo;        b >>= 1;    }    return ret;}lint loop ( lint k, lint len ){    lint ret = mod_exp(k-1, len);    if ( len & 1 ) ret = (ret + modulo - (k-1)) % modulo;    else ret = (ret + k-1) % modulo;    return ret;}//費馬小定理 a^p-1 = 1 (mod p)lint inverse ( lint n ){    return mod_exp(n,modulo-2);}struct Factor { int b, e; } f[1000];int fnum;void split ( int n ){    fnum = 0;    for ( int i = 0; i < pn && p[i] * p[i] <= n; i++ )    {        if ( n % p[i] ) continue;        f[fnum].b = p[i]; f[fnum].e = 0;        while ( n % p[i] == 0 ) { f[fnum].e++; n /= p[i]; }        fnum++;    }    if ( n > 1 )        f[fnum].b = n, f[fnum++].e = 1;}// L即是迴圈所分解的輪換的階void DFS ( int dep, int L ){    if ( dep == fnum )    {        ans = (ans + Euler(L) * loop(k-1,n/L)) % modulo;        return;    }    for ( int val = 1, i = 0; i <= f[dep].e; i++, val *= f[dep].b )        DFS ( dep + 1, L * val );}int main(){    Prime(); Euler();    while ( scanf("%I64d%I64d",&n,&k) != EOF )    {        ans = 0;        split(n);        DFS (0,1);        ans = k * ans % modulo;        ans = inverse(n) * ans % modulo;        printf("%I64d\n",ans);    }    return 0;}/* 擴充歐幾裡得求逆元lint ext_gcd ( lint a, lint b, lint& x, lint& y ){    lint ret, tmp;    if ( b == 0 )    {        x = 1, y = 0;        return a;    }    ret = ext_gcd(b, a % b, x, y);    tmp = x, x = y, y = tmp - a / b * y;    return ret;}lint inverse ( lint n ){    lint x, y;    ext_gcd (n,modulo,x,y);    x = x % modulo;    return x >= 0 ? x : x + modulo;}//用普通方法求不等價著色數lint polya ( lint n, lint k ){    lint ret = 0;    for( lint l = 1; l * l <= n; l++ )    {       if ( n % l ) continue;       ret = (ret + Euler(l) * loop(k-1, n/l)) % modulo;       if ( l * l == n ) break;       ret = (ret + Euler(n/l) * loop(k-1, l)) % modulo;    }    return ret * inverse(n) % modulo;}*/