首页 > 代码库 > 初涉矩阵快速幂
初涉矩阵快速幂
一般用来加速递推。
简单的,对于fib数列有,f0 = 1,f1 = 1,fn = fn-1 + fn-2(n >= 2)。
则对于fn有:
一般的,对于fn = A1*f(n-1) + A2*f(n-2) + .... +A(n-1)*f1,有:
又因为矩阵乘法满足结合律,所以可以用快速幂来求A^n,从而达到递推的效果。
顺便即一个小技巧:
以POJ 3233为例
#include <algorithm> #include <iostream> #include <cstring> #include <cstdlib> #include <cstdio> #include <queue> #include <cmath> #include <stack> #include <map> #pragma comment(linker, "/STACK:1024000000"); #define EPS (1e-8) #define LL long long #define ULL unsigned long long #define _LL __int64 #define INF 0x3f3f3f3f using namespace std; int Mod; const int MAXN = 61; struct Mat { LL mat[MAXN][MAXN]; int r,c; void Init(int val,int R,int C) { r = R,c = C; for(int i = 1;i <= r; ++i) for(int j = 1;j <= c; ++j) if(i != j) mat[i][j] = 0; else mat[i][j] = val; } }; Mat MatrixMult(Mat a,Mat b) { Mat p; p.Init(0,a.r,b.c); for(int i = 1;i <= a.r; ++i) { for(int j = 1;j <= b.c; ++j) { for(int k = 1;k <= b.r; ++k) { p.mat[i][j] += a.mat[i][k]*b.mat[k][j]; p.mat[i][j] %= Mod; } } } return p; } Mat QuickMult(_LL k,Mat coe) { Mat p; p.Init(1,coe.r,coe.c); while(k >= 1) { if(k&1) p = MatrixMult(p,coe); coe = MatrixMult(coe,coe); k >>= 1; } return p; } int main() { _LL n,k,m; int i,j; Mat A,B; scanf("%lld %lld %lld",&n,&k,&m); Mod = m; for(i = 1;i <= n; ++i) { for(j = 1;j <= n; ++j) scanf("%lld",&A.mat[i][j]); } for(i = 1;i <= n; ++i) { for(j = 1;j <= n; ++j) { if(i == j) A.mat[i][j+n] = 1; else A.mat[i][j+n] = 0; } } for(i = 1;i <= n; ++i) { for(j = 1;j <= n; ++j) { A.mat[i+n][j] = 0; } } for(i = 1;i <= n; ++i) { for(j = 1;j <= n; ++j) { if(i == j) A.mat[i+n][j+n] = 1; else A.mat[i+n][j+n] = 0; } } A.r = 2*n,A.c = 2*n; A = QuickMult(k+1,A); for(i = 1;i <= n; ++i) { if(A.mat[i][i+n]) A.mat[i][i+n]--; else A.mat[i][i+n] = m-1; } for(i = 1;i <= n; ++i) { for(j = 1;j <= n; ++j) { printf("%lld",A.mat[i][j+n]); if(j == n) printf("\n"); else printf(" "); } } return 0; }
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。