您的位置:首页 > 其它

图像处理之基础---矩阵求逆实现

2014-09-05 15:33 357 查看
最近做一个加密算法遇到需要计算矩阵的逆,闲着无聊,记录一下,以后免得再麻烦。

[cpp] view plaincopyprint?
#include
#include
#include

#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;
}
http://blog.csdn.net/shanshanpt/article/details/16820325
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: