a21 a22 a23 a24
a31 a32 a33 a34
a41 a42 a43 a44
假设一个二维数组中存储了以上行列式,应该怎样计算该行列式的值呀?
see http://en.wikipedia.org/wiki/Determinant
As far as I know, there are at least 3 methods:
0. Leibniz formula --- calculate all n! products and sum them up;
1. Recur on minors;
2. Gaussian elimination.
/*---------------------------------------------------------------------------
File name: determinant.c
Author: HJin (email: fish_sea_bird [at] yahoo [dot] com )
Created on: 10/29/2007 20:38:40
Environment: WinXPSP2 En Pro + VS2005 v8.0.50727.762
use the uppper-triangulation process to calculate the
determinant of a matrix.
*/
#include <stdio.h>
#include <math.h>
#define N 3
#define EPS 1e-10
void PrintMatrix(double a[][N])
{
int i, j;
for (i = 0; i < N; ++i)
{
for (j = 0; j < N; ++j)
{
printf("%8.4g ", a[i][j]);
}
printf("\n");
}
printf("\n");
}
/*
* divides a row of a by the divisor. divisor != 0.0.
*/
void DivideRow(double a[][N], int r, double divisor)
{
int j;
for (j = 0; j < N; ++j)
{
a[r][j] /= divisor;
}
}
/*
* row1 := row 1 - row2.
*/
void SubtractRow(double a[][N], int r1, int r2)
{
int j;
for (j = 0; j < N; ++j)
{
a[r1][j] -= a[r2][j];
}
}
/*
* Upper-triangulates a matrix.
* returns the determinant of the matrix.
*/
double LU(double a[][N])
{
int i, j, k;
double factor = 1.0;
double temp;
// skip row 0
for (i = 1; i < N; ++i)
{
for (j = 0; j < i; ++j)
{
if (i != j)
{
temp = a[i][j];
if (fabs(temp) > EPS)
{
DivideRow(a, i, temp);
factor *= temp;
}
}
temp = a[j][j];
if (fabs(temp) > EPS && fabs(temp - 1.0) > EPS)
{
DivideRow(a, j, temp);
factor *= temp;
}
if (fabs(a[i][j]) > EPS)
{
SubtractRow(a, i, j);
}
for (k = 0; k < N; ++k)
{
if (fabs(a[k][k]) <= EPS)
{
return 0.0;
}
}
}
}
a[N-1][N-1] *= factor;
return a[N-1][N-1];
}
double multDia(double a[][N])
{
int i;
double factor = 1.0;
for (i = 0; i < N; ++i)
{
factor *= a[i][i];
}
return factor;
}
int main()
{
double a[N][N] =
{
{8, 0, 1},
{2, 1, 3},
{5, 3, 9}
};
double factor;
PrintMatrix(a);
factor = LU(a);
PrintMatrix(a);
printf("determinant = %g.\n", factor);
return 0;
}