c语言math函数 sgn,常用矩阵计算C语言代码

参考资料：

行列式：http://zh.wikipedia.org/wiki/行列式#.E4.BB.A3.E6.95.B0.E4.BD.99.E5.AD.90.E5.BC.8F

伴随矩阵：http://zh.wikipedia.org/wiki/伴随矩阵

余因子矩阵：http://zh.wikipedia.org/wiki/余因子矩阵

逆矩阵：http://zh.wikipedia.org/wiki/逆矩阵

关于求一个矩阵的行列式，网上好多代码其实都是有问题的，我看到好多求行列式的代码都只是简单地把所有正对角线的元素乘起来再求和，然后在减去所有负对角矩阵的相应元素的乘积。这种方法在矩阵的阶大于等于4的时候是有问题的，漏掉了好多因子。正确的做法是参照行列式的定义，可以查看文章顶部的参考资料。

求矩阵行列式的正确代码如下：

//--------------------------------------------------------

//功能：求矩阵 n X n 的行列式

//入口参数：矩阵首地址 p；矩阵行数 n

//返回值：矩阵的行列式值

//--------------------------------------------------------

double determinant(double *p, int n)

{

int *list = new int[n];

for (int i = 0; i < n; i++)

list[i] = i;

double ret = det(p, n, 0, list, 0.0);

delete[] list;

return ret;

}

double det(double *p, int n, int k, int list[], double sum)

{

if(k >= n)

{

int order = inver_order(list, n);

double item = (double)sgn(order);

for (int i = 0; i < n; i++)

{

//item *= p[i][list[i]];

item *= *(p + i * n + list[i]);

}

return sum + item;

}

else

{

for(int i = k; i < n; i++)

{

swap(&list[k], &list[i]);

sum = det(p, n, k+1, list, sum);

swap(&list[k], &list[i]);

}

}

return sum;

}

void swap(int *a, int *b)

{

int m;

m = *a;

*a = *b;

*b = m;

}

//求逆序对的个数

int inver_order(int list[], int n)

{

int ret = 0;

for(int i = 1; i < n; i++)

for (int j = 0; j < i; j++)

if (list[j] > list[i])

ret++;

return ret;

}

int sgn(int order)

{

return order % 2 ? -1 : 1;

}

当然还可以用LU分解法来求，在矩阵的阶比较大时，用高斯消元法或者LU分解法求解具有一定的优势。

由于行列式是求矩阵的代数余子式的基础，代数余子式又是求矩阵的伴随矩阵的基础，求出伴随矩阵之后才可以求矩阵的逆矩阵。A矩阵的逆矩阵等于A矩阵的伴随矩阵除以A矩阵的行列式。

求矩阵某个元素的代数余子式的代码如下：

//----------------------------------------------------

//功能：求k×k矩阵中元素A(mn)的代数余子式

//入口参数：k×k矩阵首地址；元素A的下标m,n; 矩阵行数 k

//返回值： k×k矩阵中元素A(mn)的代数余子式

//----------------------------------------------------

double algebraic_cofactor(double *p, int m, int n, int k)

{

int len = (k - 1) * (k - 1);

double *cofactor = new double[len];

int count = 0;

int raw_len = k * k;

for (int i = 0; i < raw_len; i++)

if (i / k != m && i % k != n)

*(cofactor + count++) = *(p + i);

double ret = determinant(cofactor, k - 1);

if ((m + n) % 2)

ret = -ret;

delete[] cofactor;

return ret;

}

求伴随矩阵的代码如下：

//----------------------------------------------------

//功能：求k×k矩阵的伴随矩阵

//入口参数：m是k×k矩阵首地址；矩阵行数 k；输出参数adj是伴随矩阵的入口地址

//返回值： 无

//----------------------------------------------------

void adjoint_m(double *m, double *adj, int k)

{

int len = k * k;

int count = 0;

for (int i = 0; i < len; i++)

{

*(adj + count++) = algebraic_cofactor(m, i % k, i / k, k);

}

}

求逆矩阵的代码如下：

//----------------------------------------------------

//功能：求k×k矩阵的逆矩阵

//入口参数：m是k×k矩阵首地址；矩阵行数 k；输出参数inv是逆矩阵的入口地址

//返回值： 无

//----------------------------------------------------

void inverse_matrix(double *raw, double *inv, int k)

{

double det = determinant(raw, k); //求行列式

if (det == 0)

{

cout << "矩阵不可逆" << endl;

return;

}

adjoint_m(raw, inv, k); //求伴随矩阵

int len = k * k;

for (int i = 0; i < len; i++)

*(inv + i) /= det;

}

两个矩阵相乘的代码如下：

//--------------------------------------------------------

//功能：求矩阵a和b的相乘结果

//入口参数：矩阵首地址 a和b；矩阵a行数ra和列数rc；矩阵b的行数rb和列数cb

//返回值：矩阵a和b的相乘结果

//--------------------------------------------------------

double* m_multiply(double *a, double *b, double *c, int ra, int ca, int rb, int cb)

{

if (ca != rb)

{

cout << "矩阵不可乘" << endl;

return NULL;

}

double *ret = c;

if (NULL == ret)

{

ret = new double[ra * cb];

}

for (int i = 0; i < ra; i++)

for (int j = 0; j < cb; j++)

{

//double sum = a[i][0] * b[0][j];

double sum = *(a + i * ca) * (*(b + j));

for (int k = 1; k < ca; k++)

//sum += a[i][k] * b[k][j];

sum += *(a + i*ca + k) * (*(b + k*cb + j));

//c[i][j] = sum;

*(ret + i*cb + j) = sum;

}

return ret;

}

测试程序代码如下：

#include

#include "matrix.h"

using namespace std;

void main()

{

//double a[] = { 2, 6, 3, 1, 0, 2, 5, 8, 4 };

double b[] = { 1, 4, 7, 3, 0, 5, -1, 9, 11 };

double a[] = { -3, 2, -5, -1, 0, -2, 3, -4, 1 };

double det = determinant(a, 3);

cout << det << endl;

det = algebraic_cofactor(a, 1, 2, 3);

cout << det << endl;

double *adj = new double[9];

adjoint_m(a, adj, 3);

for (int i = 0; i < 9; i++)

cout << adj[i] << " ";

cout << endl;

double *c = m_multiply(a, b, NULL, 3, 3, 3, 3);

for (int i = 0; i < 9; i++)

cout << c[i] << " ";

cout << endl;

}