poj Matrix Power Series

最新推荐文章于 2020-12-19 11:24:01 发布

clover_hxy

最新推荐文章于 2020-12-19 11:24:01 发布

阅读量415

点赞数

分类专栏：矩阵

本文链接：https://blog.youkuaiyun.com/clover_hxy/article/details/50697916

版权

矩阵专栏收录该内容

27 篇文章

订阅专栏

Matrix Power Series

Time Limit: 3000MS		Memory Limit: 131072K
Total Submissions: 18956		Accepted: 8012

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A² + A³ + … + A^k.

Input

The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 10⁹) and m (m < 10⁴). 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*n的矩阵，求矩阵S = A + A² + A³ + … + A^k

^{输出对m取模的结果}

^{k的范围是10^9，直接计算时间复杂度十分不合理，所以我们可以通过递归实现时间上的缩短}

当k为偶数时，s(k)=(1+A^(k/2))*(A+A^2+…+A^(k/2))

当k为奇数时，s(k)=A+(A+A^(k/2+1))*(A+A^2+…+A^(k/2)

剩下的问题就是一些矩阵的相关运算啦

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
int n,k,p;
struct data
{
 int m[41][41];
};data a,e;//e是单位矩阵，单位矩阵有一个很好的性质A*E=E*A=A
void init(data &x,data y)
{
  for (int i=1;i<=n;i++)
   for (int j=1;j<=n;j++)
    x.m[i][j]=y.m[i][j];
} 
data add(data x,data y)//矩阵加法
{
	data c;
	for (int i=1;i<=n;i++)
	 for (int j=1;j<=n;j++)
	  c.m[i][j]=0,c.m[i][j]=(x.m[i][j]+y.m[i][j])%p;
	return c;
}
data mul(data x,data y)//矩阵乘法
{
  data temp;
  for (int i=1;i<=n;i++)
   for (int j=1;j<=n;j++)
   	{
   	  temp.m[i][j]=0;
	  for (int k=1;k<=n;k++)
	   temp.m[i][j]=(temp.m[i][j]%p+(x.m[i][k]*y.m[k][j])%p)%p;	 
   	}
  return temp;
} 
data quickpow(int x)//矩阵快速幂
{
  data sum;
  memset(sum.m,0,sizeof(sum.m));
  for (int i=1;i<=n;i++)
     sum.m[i][i]=1;
  data num;
  init(num,a);
  while (x>0)
   {
   	 if (x&1)
   	  sum=mul(sum,num);
   	 x=x>>1;
   	 num=mul(num,num);
   }
  return sum;
}
data solve(int k)
{
   if (k==1) return a;
   data b=quickpow((k+1)/2);
   data c=solve(k/2);
   if (k%2)
    return add(a,mul(add(a,b),c));
   else
    return mul(add(e,b),c); 
}
int main()
{
  scanf("%d%d%d",&n,&k,&p);
  for (int i=1;i<=n;i++)
   for (int j=1;j<=n;j++)
    scanf("%d",&a.m[i][j]);
  data ans;
  memset(e.m,0,sizeof(e.m));
  for (int i=1;i<=n;i++)
   e.m[i][i]=1;
  ans=solve(k);
  for (int i=1;i<=n;i++)
  {
   for (int j=1;j<=n;j++)
    printf("%d ",ans.m[i][j]);
   printf("\n");
  }
}