您的位置:首页 > 其它

【矩阵快速幂】POJ 3233 && NYOJ 298 Matrix Power Series

2017-04-21 15:45 429 查看
Matrix Power Series

Time Limit: 3000MS Memory Limit: 131072K
Total Submissions: 22417 Accepted: 9392
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]


Home Page   

Go
Back  

To top

转载自: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

;
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]);
4000
}
printf("\n");
}
}
return 0;
}


方法一代码:

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const int N = 32;
struct Matrix {
int mat

;
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;
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: