C实现计算任意阶行列式

#include <stdio.h>
#include <malloc.h>

typedef double DataType;
DataType CalcDeterminant(const DataType * d, int order); /* 计算行列式 */
DataType CalcDeterminantA(const DataType * d, int order, int i, int j)

/* 计算行列式的代数余子式 */
{
 int a = 0;
 int b = 0;
 int c = 0;
 DataType * d1 = (DataType *)malloc(sizeof(DataType) * (order - 1) * (order - 1));
 for(a = 0; a != order; ++a)
  if(a != i)
   for(b = 0; b != order; ++b)
    if(b != j)
     d1[c++] = d[a * order + b];
 --order;
 {
  DataType result = ((i + j) & 0X1) ? -CalcDeterminant(d1, order) : CalcDeterminant(d1, order);
  free(d1);
  return result;
 }
}
DataType CalcDeterminant(const DataType * d, int order)
{
 switch(order)
 {
 case 1: return d[0];
 case 2: return d[0] * d[3] - d[1] * d[2];
 case 3:
  {
   return d[0] * d[4] * d[8]
   + d[1] * d[5] * d[6]
   + d[2] * d[3] * d[7]
   - d[0] * d[5] * d[7]
   - d[1] * d[3] * d[8]
   - d[2] * d[4] * d[6];
  }
 default:
  {
   DataType result = 0;
   int i = 0;
   for(i = 0; i != order; ++i)
    result += d[i] * CalcDeterminantA(d, order, 0, i);
   return result;
  }
 }
}

int main(int argc, char *argv[]) /* 测试代码 */
{
 {
  double deter[] =
  {
    3, -2, 2, 1,
  };
  printf("%lf/n", CalcDeterminant(deter, 2)); /* 7 */
 }
 {
  double deter[] =
  {
   1, 2, -4,
   -2, 2, 1,
   -3, 4, -2,
  };
  printf("%lf/n", CalcDeterminant(deter, 3)); /* -14 */
 }
 {
  double deter0[] =
  {
   3,  1, -1, 2,
   -5,  1, 3, -4,
    2,  0, 1, -1,
    1, -5, 3, -3,
  };
  double deter1[] =
  {
   3, 1, 1, 1,
   1, 3, 1, 1,
   1, 1, 3, 1,
   1, 1, 1, 3,
  };
  double deter2[] =
  {
   2, 1, -5, 1,
   1, -3, 0, -6,
   0, 2, -1, 2,
   1, 4, -7, 6,
  };
  printf("%lf/n", CalcDeterminant(deter0, 4)); /* 40 */
  printf("%lf/n", CalcDeterminant(deter1, 4)); /* 48 */
  printf("%lf/n", CalcDeterminant(deter2, 4)); /* 27 */
 }
 {
  double deter0[] =
  {
   -2, 3, 0, 0, 0,
   0, -2, 3, 0, 0,
   0, 0, -2, 3, 0,
   0, 0, 0, -2, 3,
   3, 0, 0, 0, -2,
  };
  double deter1[] =
  {
   0, 3, 0, 0, 0,
   0, 0, 3, 0, 0,
   0, 0, 0, 3, 0,
   0, 0, 0, 0, 3,
   3, 0, 0, 0, 0,
  };
  printf("%lf/n", CalcDeterminant(deter0, 5)); /* 211 */
  printf("%lf/n", CalcDeterminant(deter1, 5)); /* 243 */
 }
 return 0;
}
这个程序的算法由两个函数（CalcDeterminant和CalcDeterminantA）相互递归调用实现。

用到了数学定理：

行列式等于它的任一行（列）的各元素与其对应的代数余子式乘积之和，即

D = a_i1A_i1 + a_i2A_i2 + ... + a_inA_in (i = 1, 2, ..., n),

或

D = a_1jA_1j + a_2jA_2j + ... + a_njA_nj (j = 1, 2, ..., n).

读者现在一定对此程序的算法了如指掌了吧。

呵呵，让我提一下代数余子式的定义吧。

在n阶行列式中，把(i, j)元a_ij所在的第i行和第j列划去后，留下来的n-1阶行列式叫做(i, j)元a_ij的余子式，记做M_ij；记

A_ij = (-1)^i+jM_ij,

A_ij叫做(i, j)元a_ij的代数余子式。

参考资料：

1.同济大学应用数学系.线性代数 第四版.北京：高等教育出版社，2003