题意:用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次幂为:
则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不是问题了。
#include<cmath>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<cstdio>
using namespace std;
#define lint __int64
const 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;
}
*/