当前位置: 代码迷 >> 综合 >> poj3233Matrix Power Series (等比矩阵求和)
  详细解决方案

poj3233Matrix Power Series (等比矩阵求和)

热度:92   发布时间:2023-09-30 17:25:51.0

对于整数的等比数列求和,有:

设S[n]=a+a^2+a^3+...a^n

(n==0)S[n]=1

(n为偶数)

S[n]=( a^(n/2)+1 )*S[n/2]

(n为奇数)

S[n]=( a^(n/2+1)+1 )*S[n/2]+( a^(n/2+1) )

类比上式,对于等比矩阵:

设S[n]=A+A^2+A^3+...A^n,I为单位矩阵

(n==0)S[n]=I

n为偶数和奇数时,计算公式与上面的式子相同



poj3233: 

题意: 有n*n的矩阵A,求S=(A+A^2+A^3...A^k)%m

套公式。

代码:

#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<cstdio>
#include<cmath>
#include<cstring>
#define N 33
using namespace std;
inline int read()
{int x=0,f=1;char ch=getchar();while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}return x*f;
}
int n,k,mod,I[N][N];
int a[N][N],ret[N][N];
void mul(int a[N][N],int b[N][N],int c[N][N])
{int tmp[N][N];for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)tmp[i][j]=0;for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)for(int k=1;k<=n;k++)tmp[i][j]+=a[i][k]*b[k][j]%mod,tmp[i][j]%=mod;for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)c[i][j]=(tmp[i][j]+mod)%mod;
}
void work(int a[N][N],int num,int c[N][N])
{int s[N][N];for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)s[i][j]=a[i][j];for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)c[i][j]=0;for(int i=1;i<=n;i++)c[i][i]=1;while(num){if(num&1)mul(c,s,c);mul(s,s,s);num>>=1;}
}
void cal(int a[N][N],int num,int c[N][N])
{int t[N][N],now[N][N],ret[N][N];if(num==1){for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)c[i][j]=a[i][j];return;}else if(num==0){for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)c[i][j]=a[i][j];return;}if(num&1){work(a,num/2+1,now);cal(a,num/2,ret);work(a,num/2+1,t);for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)t[i][j]+=I[i][j],t[i][j]%=mod;mul(ret,t,ret);for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)ret[i][j]+=now[i][j],ret[i][j]%=mod;for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)c[i][j]=(ret[i][j]+mod)%mod;}else{cal(a,num/2,ret);work(a,num/2,now);for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)now[i][j]+=I[i][j],now[i][j]%=mod;mul(now,ret,ret);for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)c[i][j]=(ret[i][j]+mod)%mod;}
}
/*
2 4 100
1 1
2 1
*/
void init()
{for(int i=1;i<=n;i++)I[i][i]=1;
}
int main()
{freopen("in.txt","r",stdin);freopen("my.txt","w",stdout);while(~scanf("%d",&n)){memset(a,0,sizeof(a));memset(ret,0,sizeof(ret));k=read(),mod=read();for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)a[i][j]=read();init();cal(a,k,ret);for(int i=1;i<=n;i++){for(int j=1;j<=n;j++)printf("%d ",ret[i][j]);printf("\n");}}
}


  相关解决方案