C - Matrix Power Series POJ - 3233

本文介绍了一种使用矩阵快速幂算法解决特定数学问题的方法。通过构建特定形式的矩阵,并利用快速幂的思想进行高效运算,可以求解如求矩阵的幂次等复杂计算问题。文中提供了详细的C语言实现代码。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

不多说贴代码

牛顿迭代

#include<stdio.h>
#include<string.h>
int n,m,k;

struct matrix{  
    int a[33][33];  
}a,b,c;
//s13=s6+a7+a7*s6;
//s12=s6+a6*s6;
//A+A^2+A^3+A^4+A^5+A^6=A+A^2+A^3 +A^3*(A+A^2+A^3)
//A+A^2+A^3+A^4=  A+A^2 + A^2*(A+A^2)  //偶 
//A+A^2+A^3=  A+  A^2 +A^2*(A)         //积 
//A+A^2 = A+A*A
//A


matrix cheng(matrix x,matrix y){  //x*y;
	memset(c.a,0,sizeof(c));
	int i,k,j;
	
	for(i=0;i<n;i++)
	 for(j=0;j<n;j++)
	  for(k=0;k<n;k++){
	   c.a[i][j]+=(x.a[i][k]*y.a[k][j]);	
	   if(c.a[i][j]>=m) c.a[i][j]%=m;
	  }
	  
	   
	 return c;
	 
}

matrix jia(matrix x,matrix y){  //x+y;
  	memset(c.a,0,sizeof(c));  
	int i,j;
	for(i=0;i<n;i++)
	 for(j=0;j<n;j++){
	  c.a[i][j]=(x.a[i][j] + y.a[i][j]);
	  if(c.a[i][j]>=m) c.a[i][j]%=m;	
	 }
	 
	 return c;
}


//s13=s6+a7+a7*s6;
//s12=s6+a6*s6;
//A+A^2+A^3+A^4+A^5+A^6+A^7=s3+a4+a4*s3
//A+A^2+A^3+A^4+A^5+A^6=A+A^2 +A^3 +A^3*(A+A^2+A^3)
//A+A^2+A^3+A^4+A^5=s2+a3+a3*s2//奇 
//A+A^2+A^3+A^4=  A+A^2 + A^2*(A+A^2)  //偶 
//A+A^2+A^3=  A+  A^2 +A^2*(A)         //积 
//A+A^2 = A+A*A
//A

matrix ans(matrix e, int x){
	matrix aa,bb,cc;
	int i,j; 
	if(x==1){
		return a; 
	}
	if(x&1){//奇数 
		aa=ans(e,x>>1);
		cc=cheng(a,e);
		bb=jia(aa,jia(cc,cheng(cc,aa)));
		a=cheng(cc,a);
	}
	else{//偶数 
	   aa=ans(e,x>>1);
	   bb=jia(aa,cheng(a,aa));
	   
	   a=cheng(a,a);

	}
/*	
	printf("\n\nk=%d\n",x);
	 for(i=0;i<n;i++)
	 for(j=0;j<n;j++)
	  printf("%d%c",bb.a[i][j]," \n"[j==n-1]);

*/	
	return bb;
}

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

	  		 	
		 b=ans(a,k);
 
	for(i=0;i<n;i++)
	 for(j=0;j<n;j++)
	  printf("%d%c",b.a[i][j]," \n"[j==n-1]);
	
		
	}

	return 0;
} 
矩阵套矩阵

#include<stdio.h>
#include<string.h>
int n,m,k;
struct matrix{  
    int a[33][33];  
}a,b,c;

struct node{
    matrix a[2][2];
};

/* 
    s[k]=A+A^2+...+A^k 
    矩阵: 
        |s[0] A|*|I O|^k=|s[k] A^(k+1)|  (I为单位矩阵,O为0矩阵,A为输入的矩阵) 
        |O    O| |I A|   |O    O      |   
*/  

matrix cheng(matrix x, matrix y ){
	matrix c;
	memset(c.a,0,sizeof(c.a));
	int i,j,k;
	for(i=0;i<n;i++)
	 for(j=0;j<n;j++)
	  for(k=0;k<n;k++){
	  	c.a[i][j]+=(x.a[i][k] * y.a[k][j]);
	   if(c.a[i][j]>=m) c.a[i][j]%=m; 
	  }
	  return c;
}

matrix jia(matrix x,matrix y){
	matrix c;
	memset(c.a,0,sizeof(c.a));
	int i,j;
	for(i=0;i<n;i++)
	 for(j=0;j<n;j++){
	  	c.a[i][j]=(x.a[i][j] + y.a[i][j]);
	   if(c.a[i][j]>=m) c.a[i][j]%=m; 
	  }
	return c;	
}


node  cheng_matrix(node x,node y){
	node c;
	
	int i,j,k;
	for(i=0;i<2;i++)
	 for(j=0;j<2;j++)
	 memset(c.a[i][j].a,0,sizeof(c.a[i][j].a));
	 
	 for(i=0;i<2;i++)
	 for(j=0;j<2;j++)
	 for(k=0;k<2;k++){
	 	c.a[i][j]=jia( c.a[i][j] , cheng(x.a[i][k], y.a[k][j] ));
	 }
	 
	return c;
	
	
}



int main(){
int i,j,k1,k2;
node x,y;
	while(~scanf("%d%d%d",&n,&k,&m)){
	
		memset(a.a,0,sizeof(a.a));
		for(i=0;i<n;i++){
			b.a[i][i]=1;
		for(j=0;j<n;j++){
	  scanf("%d",&(a.a[i][j]));
	  if(a.a[i][j]>=m) a.a[i][j]%=m;		 	
		 }	
		 
		}
	/* 
    s[k]=A+A^2+...+A^k 
    矩阵: 
        |s[0] A|*|I O|^k=|s[k] A^(k+1)|  (I为单位矩阵,O为0矩阵,A为输入的矩阵) 
        |O    O| |I A|   |O    O      |   
    */  
	
		
		x.a[0][0]=c;
		x.a[1][0]=c;
		x.a[1][1]=c;
		y.a[0][1]=c;
		x.a[0][1]=a;
		y.a[1][1]=a;
		
		y.a[0][0]=b;
		  y.a[1][0]=b;
		
/*			for(i=0;i<n;i++)
	 for(j=0;j<n;j++)
	  printf("%d%c",x.a[0][0].a[i][j]," \n"[j==n-1]);
	
		 printf("\n\n\n\n\n\n\n");
*/ 
	  	while(k){
	  		if(k&1) x=cheng_matrix(x,y);
			  y=cheng_matrix(y,y);
			  k>>=1; 
			  
		/*	 
			  for(k1=0;k1<2;k1++)
			   for(k2=0;k2<2;k2++){
			   	printf("X%d%d:\n",k1,k2);
			   	 for(i=0;i<n;i++)
	 for(j=0;j<n;j++)
	  printf("%d%c",x.a[k1][k2].a[i][j]," \n"[j==n-1]);
			   	
			   }
	
			  
	          for(k1=0;k1<2;k1++)
			   for(k2=0;k2<2;k2++){
			   	printf("y%d%d:\n",k1,k2);
	 for(i=0;i<n;i++)
	 for(j=0;j<n;j++)
	  printf("%d%c",y.a[k1][k2].a[i][j]," \n"[j==n-1]);
			  */
		}
			 
		  

	for(i=0;i<n;i++)
	 for(j=0;j<n;j++)
	  printf("%d%c",x.a[0][0].a[i][j]," \n"[j==n-1]);
	
		
	}
	return 0;
	
}


评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值