矩阵的求逆

我会带着你远行 2022-09-18 10:58 172阅读 0赞

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

#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;
    }