HDU Birthday Toy 特殊限制的polya
来源:互联网 发布:羊绒背心怎么样 知乎 编辑:程序博客网 时间:2024/05/16 05:16
题意:用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_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;}*/
- HDU Birthday Toy 特殊限制的polya
- hdu 2865 Birthday Toy 及我对polya的总结
- hdu 2865 Birthday Toy (polya,好题)
- HDU 2865 Birthday Toy(Polya综合)
- HDU 2865 Birthday Toy(Polya+矩阵乘法+dp)
- hdu 2481 Birthday Toy
- Birthday Toy HDU
- HDU 2865 Birthday Toy polya 矩阵快速幂 欧拉函数
- HDU 2481 Toy(Polya综合)
- HDU 5080 Colorful Toy (polya定理)
- HDU 2865 Birthday Toy(ploya好题)
- hdu 5080 - Colorful Toy(2014 AnShan)几何+polya
- HDU 5080 Colorful Toy(polya+计算几何)
- hdu 2865 Birthday Toy (置换+DP+矩阵乘法)
- HDU 2481 Toy(08成都现场 Polya,递推,矩阵,数论……)
- HDU 2089 数位DP + 特殊数字的限制
- POJ 2888 Magic Bracelet 有限制的polya
- POJ 2888 Magic Bracelet 有限制的Polya计数
- java.net.SocketException: Connection reset
- Windows7系统下取消共享文件夹后再取消共享文件夹上的小锁图标的方法
- 指针与const限定符
- Linux服务Telnet远程登录配置
- 题目1004:Median
- HDU Birthday Toy 特殊限制的polya
- 从注册表中关闭windowns自动更新
- C++与C的一些区别(自己学习中遇到的)(未完待续)
- 黑马日记第四篇——面向对象(下)
- android开发我的新浪微博客户端-阅读微博UI篇(6.1)
- poj 3264 Charm Bracelet(螺背包)
- 位操作基础篇之位操作全面总结
- 康托展开
- extern 用法,全局变量与头文件(重复定义)