Matrix Power Series
Description Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak. Input The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order. Output Output the elements of S modulo m in the same way as A is given. Sample Input 2 2 4 0 1 1 1 Sample Output 1 2 2 3 Source
POJ Monthly--2007.06.03, Huang, Jinsong
|
[Submit] [Go Back] [Status] [Discuss]
转载自:http://blog.csdn.net/lyhvoyage/article/details/39584771
一种是两次二分+递归,另外一种方法是构造矩阵;
我选的是第二种方法,因为第一种的递归掌握不好;
题意:给出一个n*n的矩阵A,求A+A^2+A^3+……+A^k mod m的结果是多少?
方法一:两次二分
Sk = A + A2 + A3 + … + Ak
=(1+A^(k/2))*(A + A^2 + A^3 + … + A^(k/2) )+{A^k}
=(1+A^(k/2))*(S(k/2) )+ {Ak} 当k为偶数时没有{Ak}
即
k%2==0: S[k]=F[k/2] (1+A[k/2]);
k%2==1: S[k]=F[k-1]+A[k];
方法二:
构造一个新矩阵B = | A E |
| 0 E |
则 B^(k+1) = | A^(k+1) A^k+A^(k-1)+……+A^2 + A + 1 |
| 0 E |
所以可以先直接用矩阵快速幂求出B^(k+1),然后用左上角的那一部分减去单位矩阵就是最后要求的矩阵。
方法二代码:
#include<cstdio>
#include<cstring>
#include<queue>
#include<algorithm>
#include<map>
using namespace std;
typedef pair<int,int > P;
typedef long long LL;
const int maxn=100;
#define N 61
int n,k,mod;
struct Matrix
{
LL mat[N][N];
Matrix()
{
memset(mat,0,sizeof(mat));
}
};
Matrix mul(Matrix a,Matrix b)
{
Matrix res;
for(int i=0; i<N; i++)
for(int j=0; j<N; j++)
{
for(int k=0; k<N; k++)
{
res.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
res.mat[i][j]%=mod;
}
}
return res;
}
Matrix pow_matrix(Matrix a,LL n)//矩阵快速幂;
{
Matrix res;
for(int i=0; i<N; i++) //初始化为单位矩阵;
res.mat[i][i]=1;
while(n!=0)
{
if(n&1)
res=mul(res,a);
a=mul(a,a);
n>>=1;
}
return res;
}
int main()
{
Matrix A;
while(~scanf("%d%d%d",&n,&k,&mod))
{
memset(A.mat,0,sizeof(A.mat));
for(int i=0; i<n; i++)
for(int j=0; j<n; j++)
{
scanf("%d",&A.mat[i][j]);
A.mat[i][j]%=mod;
A.mat[i][i+n]=1;
}
for(int i=n; i<2*n; i++)
A.mat[i][i]=1;
Matrix ans=pow_matrix(A,k+1);
for(int i=0; i<n; i++)
{
for(int j=n; j<n*2; j++)
{
if(j>n)
printf(" ");
if(j-i==n)
printf("%d",(ans.mat[i][j]-1+mod)%mod);
else
printf("%d",ans.mat[i][j]);
}
printf("\n");
}
}
return 0;
}
方法一代码:
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
const int N = 32;
struct Matrix {
int mat[N][N];
Matrix() {
memset(mat, 0, sizeof(mat));
for(int i = 0; i < N; i++)
mat[i][i] = 1;
}
} E;
int n, k, mod;
Matrix Multi(Matrix a, Matrix b) {
Matrix res;
for(int i = 0; i < n; i++) {
for(int j = 0; j < n; j++) {
res.mat[i][j] = 0;
for(int k = 0; k < n; k++)
res.mat[i][j] = (res.mat[i][j] + a.mat[i][k] * b.mat[k][j]) % mod;
}
}
return res;
}
Matrix Add(Matrix a, Matrix b) {
Matrix res;
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
res.mat[i][j] = (a.mat[i][j] + b.mat[i][j]) % mod;
return res;
}
Matrix Pow(Matrix a, int n) {
Matrix res;
while(n) {
if(n&1) res = Multi(res, a);
a = Multi(a, a);
n >>= 1;
}
return res;
}
Matrix Get_Ans(Matrix a, int k) {
if(k == 1) return a;
if(k&1) return Add(Pow(a, k), Get_Ans(a, k-1));
if(k % 2 == 0) {
Matrix A = Get_Ans(a, k/2);
Matrix B = Pow(a, k/2);
Matrix C = Multi(A, B);
return Add(C, A);
}
}
int main() {
Matrix A;
while(~scanf("%d%d%d", &n, &k, &mod)) {
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++) {
scanf("%d", &A.mat[i][j]);
A.mat[i][j] %= mod;
}
Matrix ans = Get_Ans(A, k);
for(int i = 0; i < n; i++) {
for(int j = 0; j < n; j++) {
if(j) printf(" ");
printf("%d", ans.mat[i][j]);
}
printf("\n");
}
}
return 0;
}
http://www.matrix67.com/blog/archives/276;