您的位置：首页 > 其它

POJ 3233 Matrix Power Series 解题报告（子矩阵构造+矩阵快速幂）

2014-04-24 20:11 253 查看

Matrix Power Series

Time Limit: 3000MS		Memory Limit: 131072K
Total Submissions: 14105		Accepted: 6078

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

解题报告：求矩阵的n次方和。

如果接触过矩阵，都知道矩阵快速幂。一开始我也这么做的，二分计算前半段和后半段，即(A^1+A^2+...+A^(n/2))*(A^(n/2)+I)。题目是A了，但是耗时1600MS+。

看了一下Status，发现一堆0MS的……又想到整数的n次方和是有公式的，那么矩阵有吗？

百度了一下解题报告，发现了一个非常巧妙的方法：（引用：http://www.cnblogs.com/rainydays/archive/2011/02/21/1960189.html）

把问题转化以加速，令

B = A I

      0 I

则B^(k + 1) = A^(k + 1)      I + A + A2 + A3 + … + Ak

                            0                          I

用二分法求B^(k + 1)
根据该性质，复杂度大大减小。另外矩阵相乘时每次最后再取模，会快上很多。

代码如下：

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

typedef long long LL;
const int SIZE = 60;
int n, mod;

struct Matrix
{
int a[SIZE][SIZE];
Matrix(int t=0)
{
memset(a, 0, sizeof(a));
for(int i=0;i<SIZE;i++) a[i][i]=t;
}

Matrix operator*(const Matrix& b) const
{
Matrix c;
for(int i=0;i<n;i++) for(int j=0;j<n;j++)
{
LL sum=0;
for(int k=0;k<n;k++) sum+=a[i][k]*b.a[k][j];
c.a[i][j]=sum%mod;
}
return c;
}
};

Matrix pow(Matrix a, int b)
{
Matrix res(1);
while(b)
{
if(b&1)
res=res*a;
a=a*a;
b>>=1;
}
return res;
}

int main()
{
int k;
while(~scanf("%d%d%d", &n, &k, &mod))
{
Matrix a;
for(int i=0;i<n;i++) for(int j=0;j<n;j++)
scanf("%d", &a.a[i][j]), a.a[i][j]%=mod;

for(int i=0;i<n;i++)
a.a[i][i+n]=a.a[i+n][i+n]=1;

n<<=1;
a = pow(a, k+1);
n>>=1;

for(int i=0;i<n;i++)
if(a.a[i][i+n]==0)
a.a[i][i+n]=mod-1;
else
a.a[i][i+n]--;

for(int i=0;i<n;i++, puts("")) for(int j=0;j<n;j++)
printf("%d ", a.a[i][j+n]);
}
}

同样附上二分的代码：

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

int n, mod;

struct Matrix
{
int a[30][30];
Matrix(int t=0)
{
memset(a, 0, sizeof(a));
for(int i=0;i<30;i++) a[i][i]=t;
}

Matrix operator*(const Matrix& b) const
{
Matrix c;
for(int i=0;i<n;i++) for(int j=0;j<n;j++)
for(int k=0;k<n;k++) c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j])%mod;
return c;
}

Matrix operator+(const Matrix& b) const
{
Matrix c;
for(int i=0;i<n;i++) for(int j=0;j<n;j++)
c.a[i][j]=(a[i][j]+b.a[i][j])%mod;
return c;
}
} one(1);

Matrix pow(Matrix a, int b)
{
Matrix res(1);
while(b)
{
if(b&1)
res=res*a;
a=a*a;
b>>=1;
}
return res;
}

Matrix powSum(Matrix a, int b)
{
if(b==0)
return one;
else if(b==1)
return a;
else if(b%2==0)
return powSum(a, b/2)*(pow(a, b/2)+one);
else
return powSum(a, b/2)*(pow(a, b/2)+one)+pow(a, b);
}

int main()
{
int k;
while(~scanf("%d%d%d", &n, &k, &mod))
{
Matrix a;
for(int i=0;i<n;i++) for(int j=0;j<n;j++)
scanf("%d", &a.a[i][j]), a.a[i][j]%=mod;

a = powSum(a, k);
for(int i=0;i<n;i++) for(int j=0;j<n;j++)
printf("%d", a.a[i][j]%mod), printf(j==n-1?"\n":" ");
}
}

同样的，一个数的n次方和再模一个数也可以这样做。公式求法中需要除以一个数，而这个数在模mod系下不一定有逆元，所以用矩阵二分来做是个不错的方法。（和队友JX讨论的）。

内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理

标签： 矩阵 n次方和

相关文章推荐

新的分享

章节导航