矩阵的求逆

来源:互联网 发布:速卖通产品数据分析 编辑:程序博客网 时间:2024/05/01 18:55


最近做一个加密算法遇到需要计算矩阵的逆,闲着无聊,记录一下,以后免得再麻烦。


#include <stdio.h>#include <math.h>#include <string.h>#define MAX 20#define E 0.000000001/** * 计算矩阵src的模 */double calculate_A( double src[][MAX], int n ){int i,j,k,x,y;double tmp[MAX][MAX], t;double result = 0.0;if( n == 1 ){return src[0][0];}for( i = 0; i < n; ++i ){for( j = 0; j < n - 1; ++j ){for( k = 0; k < n - 1; ++k ){x = j + 1;y = k >= i ? k + 1 : k;tmp[j][k] = src[x][y];}}t = calculate_A( tmp, n - 1 );if( i % 2 == 0 ){result += src[0][i] * t;}else{result -= src[0][i] * t;}}return result;}/** * 计算伴随矩阵 */void calculate_A_adjoint( double src[MAX][MAX], double dst[MAX][MAX], int n ){int i, j, k, t, x, y;double tmp[MAX][MAX];if( n == 1 ){dst[0][0] = 1;return;}for( i = 0; i < n; ++i ){for( j = 0; j < n; ++j ){for( k = 0; k < n - 1; ++k ){for( t = 0; t < n - 1; ++t ){x = k >= i ? k + 1 : k ;y = t >= j ? t + 1 : t;tmp[k][t] = src[x][y];}}dst[j][i]  =  calculate_A( tmp, n - 1 );if( ( i + j ) % 2 == 1 ){dst[j][i] = -1*dst[j][i];}}}}/** * 得到逆矩阵 */int calculate_A_inverse( double src[MAX][MAX], double dst[MAX][MAX], int n ){double A = calculate_A( src, n );double tmp[MAX][MAX];int i, j;if ( fabs( A - 0 ) <= E ){printf("不可能有逆矩阵!\n");return 0;}calculate_A_adjoint( src, tmp, n );  for( i = 0; i < n; ++i )  {  for( j = 0; j < n; ++j )  {  dst[i][j] = (double)( tmp[i][j] / A );}  }return 1;}/** * 输出矩形查看 */void print_M( double M[][MAX], int n ){int i, j;for ( i = 0; i < n; ++i ){for ( j = 0; j < n; ++j ){printf("%lf ", M[i][j]);}printf("\n");}}/** * main */int main(){/** * 测试矩阵 */double test[MAX][MAX], dst[MAX][MAX];int n = 3;int is_exist;/** * 构造最简单的: * 1, 0, 0, *  0, 2, 0, *  0, 0, 5 */memset( test, 0, sizeof( test ) );test[0][0] = 1.0;test[1][1] = 2.0;test[2][2] = 5.0;is_exist = calculate_A_inverse( test, dst, n );if ( is_exist ){print_M(dst, n);}else{printf("不存在!\n");}return 0;}


原创粉丝点击