首页 > 代码库 > hdu 4549 M斐波那契数列(矩阵快速幂,快速幂降幂)
hdu 4549 M斐波那契数列(矩阵快速幂,快速幂降幂)
http://acm.hdu.edu.cn/showproblem.php?pid=4549
f[0] = a^1*b^0%p,f[1] = a^0*b^1%p,f[2] = a^1*b^1%p.....f[n] = a^fib[n-1] * b^fib[n-2]%p。
这里p是质数,且a,p互素,那么我们求a^b%p,当b很大时要对b降幂。
因为a,p互素,那么由费马小定理知a^(p-1)%p = 1。令b = k*(p-1) + b‘,a^b%p = a^(k*(p-1)+b‘)%p = a^(k*(p-1))%p * a^b‘%p
= a^b‘%p。其中p‘ = b%(p-1)。
那么上题的fib[n-1]和fib[n-2]可以经过矩阵快速幂对(p-1)取模,实现降幂以后再用快速幂求得。
其实a^b%c较一般的形式是a^(b%phi[c]+phi[c])%c (b>=phi[c]),在c是素数或a与c互素的时候可以简化为a^(b%(c-1))%c。
纳闷了一中午,结果把fib数列弄错了。f[0] = 0, f[1] = 1,f[n] = f[n-1]+f[n-2](n >= 2)。第0项是0,sad....
#include <stdio.h> #include <iostream> #include <map> #include <set> #include <list> #include <stack> #include <vector> #include <math.h> #include <string.h> #include <queue> #include <string> #include <stdlib.h> #include <algorithm> #define LL long long #define _LL __int64 #define eps 1e-12 #define PI acos(-1.0) #define C 240 #define S 20 using namespace std; const int maxn = 110; const int mod = 1000000007; struct matrix { LL mat[3][3]; void init() { memset(mat,0,sizeof(mat)); for(int i = 0; i < 2; i++) mat[i][i] = 1; } }m; LL a,b,n; matrix mul_matrix(matrix x, matrix y) { matrix ans; memset(ans.mat,0,sizeof(ans.mat)); for(int i = 0; i < 2; i++) { for(int k = 0; k < 2; k++) { if(x.mat[i][k] == 0) continue; for(int j = 0; j < 2; j++) ans.mat[i][j] = (ans.mat[i][j] + x.mat[i][k]*y.mat[k][j])%(mod-1); } } return ans; } matrix pow_matrix(matrix m, LL n) { matrix ans; ans.init(); while(n) { if(n&1) ans = mul_matrix(ans,m); m = mul_matrix(m,m); n >>= 1; } return ans; } LL pow_mod(LL a, LL b) { LL res = 1; a = a%mod; while(b) { if(b&1) res = res*a%mod; a = a*a%mod; b >>= 1; } return res; } int main() { while(~scanf("%lld %lld %lld",&a,&b,&n)) { m.mat[0][0] = m.mat[0][1] = m.mat[1][0] = 1; m.mat[1][1] = 0; if(n == 0) { printf("%lld\n",a%mod); continue; } if(n == 1) { printf("%lld\n",b%mod); continue; } else if(n == 2) { printf("%lld\n",a*b%mod); continue; } matrix ans = pow_matrix(m,n); LL res = pow_mod(a,ans.mat[1][1]) * pow_mod(b,ans.mat[0][1]) %mod; printf("%lld\n",res); } return 0; }
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。