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 = ai1Ai1 + ai2Ai2 + ... + ainAin (i = 1, 2, ..., n),

或

D = a1jA1j + a2jA2j + ... + anjAnj (j = 1, 2, ..., n).

讀者現在一定對此程式的演算法了如指掌了吧。

呵呵，讓我提一下代數餘子式的定義吧。

在n階行列式中，把(i, j)元aij所在的第i行和第j列划去後，留下來的n-1階行列式叫做(i, j)元aij的餘子式，記做Mij；記

Aij = (-1)i+jMij,

Aij叫做(i, j)元aij的代數餘子式。

參考資料：

1.同濟大學應用數學系.線性代數 第四版.北京：高等教育出版社，2003