|
楼主 |
发表于 2014-11-8 16:01
|
显示全部楼层
#include <stdio.h>
signed int signmir(int *j,int n){
int i[n];
signed int signindex=1;
if(j[1]<j[0])signindex*=-1;
if (n==2)return signindex;
for (i[2]=0;i[2]<3;i[2]++)
if(j[2]<j[i[2]])signindex*=-1;
if (n==3) return signindex;
for (i[3]=0;i[3]<4;i[3]++)
if(j[3]<j[i[3]])signindex*=-1;
if (n==4) return signindex;
for (i[4]=0;i[4]<5;i[4]++)
if(j[4]<j[i[4]])signindex*=-1;
if (n==5) return signindex;
for (i[5]=0;i[5]<6;i[5]++)
if(j[5]<j[i[5]])signindex*=-1;
if (n==6) return signindex;
for (i[6]=0;i[6]<7;i[6]++)
if(j[6]<j[i[6]])signindex*=-1;
if (n==7) return signindex;
for (i[7]=0;i[7]<8;i[7]++)
if(j[7]<j[i[7]])signindex*=-1;
if (n==8) return signindex;
for (i[8]=0;i[8]<9;i[8]++)
if(j[8]<j[i[8]])signindex*=-1;
if (n==9) return signindex;
for (i[9]=0;i[9]<10;i[9]++)
if(j[9]<j[i[9]])signindex*=-1;
if (n==10) return signindex;
}
double compu( int n ,double a[10][10]){
int j[n];//计数第二个下标
double mir=0;
for(j[0]=0;j[0]<n;j[0]++){
if(n==1){mir=a[0][0];//一维
break;
}
for(j[1]=0;j[1]<n;j[1]++){
if(j[1]==j[0])continue;
if(n==2){//二维
mir+=signmir(j,n)*a[0][j[0]]*a[1][j[1]];
break;
}
for(j[2]=0;j[2]<n;j[2]++){
if((j[2]==j[0])||(j[2]==j[1]))continue;
if(n==3){//三维
// printf("in 3,a[0][%d]=%lf,a[1][%d]=%lf,a[2][%d]=%lf,\n",j[0],a[0][j[0]],j[1],a[1][j[1]],j[2],\
// a[2][j[2]]);
// printf("signmir=%d\n",signmir(j,n));
// printf("本次乘积为%lf\n",signmir(j,n)*a[0][j[0]]*a[1][j[1]]*a[2][j[2]]);
mir+=signmir(j,n)*a[0][j[0]]*a[1][j[1]]*a[2][j[2]];
// printf("mir=%lf\n",mir);
break;
}
for(j[3]=0;j[3]<n;j[3]++){
// printf("in 4,j[0]=%d,j[1]=%d,j[2]=%d\n",j[0],j[1],j[2]);
if((j[3]==j[0])||(j[3]==j[1])||(j[3]==j[2]))continue;
if(n==4){//四维
//printf("n==4");
mir+=signmir(j,n)*a[0][j[0]]*a[1][j[1]]*a[2][j[2]]*a[3][j[3]];
continue;
}
for(j[4]=0;j[4]<n;j[4]++){
// printf("in 5,j[0]=%d,j[1]=%d,j[2]=%d\n",j[0],j[1],j[2]);
if((j[4]==j[0])||(j[4]==j[1])||(j[4]==j[2])||(j[4]==j[3]))continue;
if(n==5){//五维
// printf("n==5");
// printf("in ,a[0][%d]=%lf,a[1][%d]=%lf,a[2][%d]=%lf,a[3][%d]=%lf,a[4][%d]=%lf\n",\
// j[0],a[0][j[0]],j[1],a[1][j[1]],j[2],a[2][j[2]],j[3],a[3][j[3]],j[4],a[4][j[4]]);
// a[2][j[2]]);
// printf("signmir=%d\n",signmir(j,n));
//printf("本次乘积为%lf\n",signmir(j,n)*a[0][j[0]]*a[1][j[1]]*a[2][j[2]]);
mir+=signmir(j,n)*a[0][j[0]]*a[1][j[1]]*a[2][j[2]]*a[3][j[3]]*a[4][j[4]];
continue;
}
for(j[5]=0;j[5]<n;j[5]++){
// printf("in 6,j[0]=%d,j[1]=%d,j[2]=%d\n",j[0],j[1],j[2]);
if((j[5]==j[0])||(j[5]==j[1])||(j[5]==j[2])||(j[5]==j[3])||(j[5]==j[4]))continue;
if(n==6){//六维
// printf("n==6");
mir+=signmir(j,n)*a[0][j[0]]*a[1][j[1]]*a[2][j[2]]*a[3][j[3]]*a[4][j[4]]*a[5][j[5]];
continue;
}
for(j[6]=0;j[6]<n;j[6]++){
// printf("in 7,j[0]=%d,j[1]=%d,j[2]=%d\n",j[0],j[1],j[2]);
if((j[6]==j[0])||(j[6]==j[1])||(j[6]==j[2])||(j[6]==j[3])||(j[6]==j[4])||(j[6]==j[5]))continue;
if(n==7){//七维
// printf("n==7");
mir+=signmir(j,n)*a[0][j[0]]*a[1][j[1]]*a[2][j[2]]*a[3][j[3]]*a[4][j[4]]*a[5][j[5]]*a[6][j[6]];
continue;
}
for(j[7]=0;j[7]<n;j[7]++){
// printf("in 8,j[0]=%d,j[1]=%d,j[2]=%d\n",j[0],j[1],j[2]);
if((j[7]==j[0])||(j[7]==j[1])||(j[7]==j[2])||(j[7]==j[3])||(j[7]==j[4])\
||(j[7]==j[5])||(j[7]==j[6]))continue;
if(n==8){//八维
// printf("n==7");
mir+=signmir(j,n)*a[0][j[0]]*a[1][j[1]]*a[2][j[2]]*a[3][j[3]]*a[4][j[4]]*a[5][j[5]]\
*a[6][j[6]]*a[7][j[7]];
continue;
}
for(j[8]=0;j[8]<n;j[8]++){
// printf("in 9,j[0]=%d,j[1]=%d,j[2]=%d\n",j[0],j[1],j[2]);
if((j[8]==j[0])||(j[8]==j[1])||(j[8]==j[2])||(j[8]==j[3])||(j[8]==j[4])\
||(j[8]==j[5])||(j[8]==j[6])||(j[8]==j[7]))continue;
if(n==9){//九维
// printf("n==9");
mir+=signmir(j,n)*a[0][j[0]]*a[1][j[1]]*a[2][j[2]]*a[3][j[3]]*a[4][j[4]]*a[5][j[5]]\
*a[6][j[6]]*a[7][j[7]]*a[8][j[8]];
continue;
}
for(j[9]=0;j[9]<n;j[9]++){
// printf("in 10,j[0]=%d,j[1]=%d,j[2]=%d\n",j[0],j[1],j[2]);
if((j[9]==j[0])||(j[9]==j[1])||(j[9]==j[2])||(j[9]==j[3])||(j[9]==j[4])\
||(j[9]==j[5])||(j[9]==j[6])||(j[9]==j[7])||(j[9]==j[8]))continue;
if(n==10){//十维
//printf("n==10");
mir+=signmir(j,n)*a[0][j[0]]*a[1][j[1]]*a[2][j[2]]*a[3][j[3]]*a[4][j[4]]*a[5][j[5]]\
*a[6][j[6]]*a[7][j[7]]*a[8][j[8]]*a[9][j[9]];
}
}
}
}
}
}
}
}
}
}
}
//ch=getchar();
return mir;
}
void yizhishi(double a[10][10], double b[10][10], int r, int s, int n){
//复制r,s之前的行和列 |||直接复制
int i,j;
for(i=0;i<r+1;i++){
for(j=0;j<s+1;j++){
b[i][j]=a[i][j];
}
}
//复制r之后,s之前的行和列 ||| r-1,S直接复制
for(i=r+1;i<n;i++){
for(j=0;j<s+1;j++){
b[i-1][j]=a[i][j];
}
}
//复制r之前,s之后的行和列 ||| r直接复制,S-1
for(i=0;i<r+1;i++){
for(j=s+1;j<n;j++){
b[i][j-1]=a[i][j];
}
}
//复制r之后,s之后的行和列 ||| r-1,S-1直接复制
for(i=r+1;i<n;i++){
for(j=s+1;j<n;j++){
b[i-1][j-1]=a[i][j];
}
}
}
int main(void){
int n=0;int r,s;
char ch;
printf("程序发布者:遐逸。2014-11-07 12:20\n如需超过10阶以上的逆矩阵函数,请联系QQ:103182681。\n\n");
printf("输入一个非'#'的字符(输入 '#' 将结束程序。)\n");
while ((ch=getchar())!='#'){
printf("输入一个数字在1到10之间\nn=");
scanf("%d",&n);
ch=getchar();
while ((n>10) || (n<1))
{
printf("因10时,计算量已经是以3628800为基本单位,运算共需至少11单位\n超过10之外,数据运算量十分大。\n本程序要求数字须在1到10之间,请重新输入。\nn=");
scanf("%d",&n);
ch=getchar();
}
double a[10][10];//原矩阵
int j[n];//计数第二个下标
double mir;
printf("以前面指定的%d为一行,以逗号','分割\n",n);
for (r=0;r<n;r++){
scanf("%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf",&a[r][0],&a[r][1],&a[r][2]\
,&a[r][3],&a[r][4],&a[r][5],&a[r][6],&a[r][7],&a[r][8],&a[r][9]);
}
printf("原矩阵\n");
for (r=0;r<n;r++){
for (s=0;s<n;s++){
printf("%10lf",a[r][s]);
}
printf("\n");
}
mir=compu(n,a);//返回原矩阵的行列式的值
printf("原矩阵对应行列式的值=%lf\n",mir);
if(mir==0){
printf("矩阵对应的行列式为0,无逆矩阵。\n\n");
ch=getchar();
continue;
}
//计算余子式
//1、构造出Mij行列式b[][];2、计算出Mij的值(*(-1)^(i+j));3、计算出相应的逆矩阵式 c[][]
double b[10][10],c[10][10];
int p,q;
if(n==1){
c[0][0]=1/mir;
}
else{
for (r=0;r<n;r++){
for(s=0;s<n;s++){
yizhishi(a,b,r,s,n);//构造成Mij行列式b[][];
/* printf("r=%d,s=%d\n",r,s);
for (p=0;p<n-1;p++){
for (q=0;q<n-1;q++){
printf("%10f",b[p][q]);
}
printf("\n");
}
*/
c[s][r]=compu(n-1,b)/mir;//算出逆矩阵 ;
if((r+s)%2==1)c[s][r]*=-1;//给出正负号;
// printf("c[%d][%d]=%lf, compu(%d-1=%d,b)=%lf,mir=%lf, r+s%2=%d\n",r,s,c[r][s],n,n-1,compu(n-1,b),mir,(r+s)%2);
}
}
}
printf("逆矩阵\n(如果等待时间较长,是正常现象,因为运算量较大)\n\n");
for (r=0;r<n;r++){
for (s=0;s<n;s++){
printf("%10f",c[r][s]);
}
printf("\n");
}
printf("\n验算矩阵,以对角线元素为1,其余元素为0表示求逆成功\n如果数值相差较大,则是精度不足导致。\n\n");
//验算逆矩阵
double d[10][10]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
for (r=0;r<n;r++){
for (s=0;s<n;s++){
for (p=0;p<n;p++){
d[r][s] +=a[r][p]* c[p][s];
//printf("a[%d][%d]=%lf, b[%d][%d]=%lf\n",i,p,a[i][p],p,j,b[p][j]);
}
printf("[%d][%d]=%lf, ",r,s,d[r][s]);
}
printf("\n");
}
printf("输入一个非'#'的字符(输入 '#' 将结束程序。)\n");
ch=getchar();
}
return 0;
} |
|