计算行列式的值

遇到一个问题，就是计算方块矩阵的行列式值。

在网上找了一些现成的，但是问题是太慢了，原因是这些现成的都是依据行列式值定义来进行计算的。

没办法自己按照书上的比较快速的计算方法写了个程序，如下。

基本原理就是变成上三角为0的来进行计算。

这里用到了double的原因是因为涉及到多次的乘除法消去一行，可能会有舍入误差，所以用double来尽量减少这个误差。

算法的来源见 https://books.google.com/books?id=1BxMryiz0_QC&pg=PA110&lpg=PA110&dq=%E8%A1%8C%E5%88%97%E5%BC%8F+%E8%AE%A1%E7%AE%97+%E5%BF%AB%E9%80%9F&source=bl&ots=JwUhB_F5zj&sig=jmLAjRdMJezkOFLs2YP0grAspkA&hl=zh-CN&sa=X&ved=0ahUKEwiTjrvXzsLMAhVHG5QKHe-eD-wQ6AEIaTAJ#v=onepage&q=%E8%A1%8C%E5%88%97%E5%BC%8F%20%E8%AE%A1%E7%AE%97%20%E5%BF%AB%E9%80%9F&f=false

代码如下：

#include "stdafx.h"#include "stdio.h"#include "stdlib.h"#include "math.h"#define g_order 5#define M(m, i, j) (m[(i) * g_order + (j)])int findNoneZero(double* matrix, int column);double det(double* matrix, int startRowColumn);int allSign = 1;float result = 1;void display(double* matrix) {printf("\n");for (int i=0; i<g_order; i++) {for (int j=0; j<g_order; j++) {printf("%f, ", M(matrix, i, j));}}}int main(){//initdouble a[g_order * g_order];for (int i=0; i<g_order * g_order; i++) {a[i] = (int)(20.0*rand()/(RAND_MAX+1.0));}display(a);double r = det(a, 0);printf("Result = %lf\n", r);getchar();} void exchange(double* matrix, int startRowColumn, int index) {for (int i=startRowColumn; i<g_order; i++) {double temp = M(matrix, i, index);M(matrix, i, index) = M(matrix, i, startRowColumn);M(matrix, i, startRowColumn) = temp;}}double det(double* matrix, int startRowColumn) {double result = 1;//处理第i行for (int i=startRowColumn; i<g_order; i++) {if (i == g_order - 1) {result = result * M(matrix, g_order - 1, g_order - 1);break;}int index = findNoneZero(matrix, i);if (index == -1) {return 0;}if (index != i) {exchange(matrix, startRowColumn, index);//交换了两列，变号allSign *= -1;}//现在M(i, i) != 0了// display(matrix);result *= M(matrix, i, i);double a = M(matrix, i, i);for (int j = i + 1; j < g_order; j++) {M(matrix, i, j) = M(matrix, i, j) / a;} //display(matrix); //printf("extract a = %f\n", a);//然后把矩阵内其它元素变换for (int j = i + 1; j<g_order; j++) {for (int k = i + 1; k<g_order; k++) {M(matrix, j, k) = M(matrix, j, k)          - M(matrix, j, i) * M(matrix, i, k);//printf(" set (%d, %d) = %f\n", j, k, M(matrix, j ,k));}}// printf("Last = \n");// display(matrix);}return allSign * result;}int findNoneZero(double* matrix, int row) {for (int i=row; i<g_order; i++) {if (M(matrix, row, i) != 0) {return i;}}return -1;}

另：在计算的过程中发现，double还是不够用，一些精度丢失了。要解决这个问题，只能用高精度乘法除法来计算。于是从网上找了个大整数算法，改进了上面的程序：

#include "stdafx.h"#include "stdio.h"#include "stdlib.h"#include "math.h"#defineMAXDIGITS1024/* maximum length bignum */ #define PLUS1/* positive sign bit */#define MINUS-1/* negative sign bit */typedef struct {    char digits[MAXDIGITS];         /* represent the number */int signbit;/* 1 if positive, -1 if negative */     int lastdigit;/* index of high-order digit */} bignum;void subtract_bignum(bignum *a, bignum *b, bignum *c);void multiply_bignum(bignum *a, bignum *b, bignum *c);void zero_justify(bignum *n);int compare_bignum(bignum *a, bignum *b);void print_bignum(bignum *n);#define g_order 4struct fractor {bignum top;bignum down;} ;typedef struct fractor fractor;#define M(m, i, j) (m[(i) * g_order + (j)])int findNoneZero(fractor* matrix, int column);bignum det(fractor* matrix, int startRowColumn);int allSign = 1;void display(fractor* matrix) {printf("\n");for (int i=0; i<g_order; i++) {for (int j=0; j<g_order; j++) {print_bignum(&M(matrix, i, j).top);printf("/");print_bignum(&M(matrix, i, j).down);printf(" ");}printf("\n");}}void copy_bignum(bignum *to, bignum* from){memcpy(to, from, sizeof(bignum));}void print_bignum(bignum *n){int i;if (n->signbit == MINUS) printf("-");for (i=n->lastdigit; i>=0; i--)printf("%c",'0'+ n->digits[i]);}float bignum_to_float(bignum *n){int i;float r = 0;// if (n->signbit == MINUS) r = -1.0;int first = 1;for (i=n->lastdigit; i>=0; i--) {if (first == 1) {r = n ->digits [i];first = 0;} else {r = 10 * r + n->digits[i];}}return r;}void int_to_bignum(int s, bignum *n){int i;/* counter */int t;/* int to work with */if (s >= 0) n->signbit = PLUS;else n->signbit = MINUS;for (i=0; i<MAXDIGITS; i++) n->digits[i] = (char) 0;n->lastdigit = -1;t = abs(s);while (t > 0) {n->lastdigit ++;n->digits[ n->lastdigit ] = (t % 10);t = t / 10;}if (s == 0) n->lastdigit = 0;}void initialize_bignum(bignum *n){int_to_bignum(0,n);}int max(int a, int b){if (a > b) return(a); else return(b);}/*c = a +-/* b;*/void add_bignum(bignum *a, bignum *b, bignum *c){int carry;/* carry digit */int i;/* counter */initialize_bignum(c);if (a->signbit == b->signbit) c->signbit = a->signbit;else {if (a->signbit == MINUS) {a->signbit = PLUS;subtract_bignum(b,a,c);a->signbit = MINUS;} else {            b->signbit = PLUS;            subtract_bignum(a,b,c);            b->signbit = MINUS;}return;}c->lastdigit = max(a->lastdigit,b->lastdigit)+1;carry = 0;for (i=0; i<=(c->lastdigit); i++) {c->digits[i] = (char) (carry+a->digits[i]+b->digits[i]) % 10;carry = (carry + a->digits[i] + b->digits[i]) / 10;}zero_justify(c);}void subtract_bignum(bignum *a, bignum *b, bignum *c){int borrow;/* has anything been borrowed? */int v;/* placeholder digit */int i;/* counter */initialize_bignum(c);if ((a->signbit == MINUS) || (b->signbit == MINUS)) {                b->signbit = -1 * b->signbit;                add_bignum(a,b,c);                b->signbit = -1 * b->signbit;return;        }if (compare_bignum(a,b) == PLUS) {subtract_bignum(b,a,c);c->signbit = MINUS;return;}        c->lastdigit = max(a->lastdigit,b->lastdigit);        borrow = 0;        for (i=0; i<=(c->lastdigit); i++) {v = (a->digits[i] - borrow - b->digits[i]);if (a->digits[i] > 0)borrow = 0;if (v < 0) {v = v + 10;borrow = 1;}                c->digits[i] = (char) v % 10;        }zero_justify(c);}int compare_bignum(bignum *a, bignum *b){int i;/* counter */if ((a->signbit == MINUS) && (b->signbit == PLUS)) return(PLUS);if ((a->signbit == PLUS) && (b->signbit == MINUS)) return(MINUS);if (b->lastdigit > a->lastdigit) return (PLUS * a->signbit);if (a->lastdigit > b->lastdigit) return (MINUS * a->signbit);for (i = a->lastdigit; i>=0; i--) {if (a->digits[i] > b->digits[i]) return(MINUS * a->signbit);if (b->digits[i] > a->digits[i]) return(PLUS * a->signbit);}return(0);}void zero_justify(bignum *n){while ((n->lastdigit > 0) && (n->digits[ n->lastdigit ] == 0))n->lastdigit --;        if ((n->lastdigit == 0) && (n->digits[0] == 0))n->signbit = PLUS;/* hack to avoid -0 */}void digit_shift(bignum *n, int d)/* multiply n by 10^d */{int i;/* counter */if ((n->lastdigit == 0) && (n->digits[0] == 0)) return;for (i=n->lastdigit; i>=0; i--)n->digits[i+d] = n->digits[i];for (i=0; i<d; i++) n->digits[i] = 0;n->lastdigit = n->lastdigit + d;}void multiply_bignum(bignum *a, bignum *b, bignum *c){bignum row;/* represent shifted row */bignum tmp;/* placeholder bignum */int i,j;/* counters */initialize_bignum(c);row = *a;for (i=0; i<=b->lastdigit; i++) {for (j=1; j<=b->digits[i]; j++) {add_bignum(c,&row,&tmp);*c = tmp;}digit_shift(&row,1);}c->signbit = a->signbit * b->signbit;zero_justify(c);}void divide_bignum(bignum *a, bignum *b, bignum *c){    bignum row;                     /* represent shifted row */    bignum tmp;                     /* placeholder bignum */int asign, bsign;/* temporary signs */    int i,j;                        /* counters */initialize_bignum(c);c->signbit = a->signbit * b->signbit;asign = a->signbit;bsign = b->signbit;a->signbit = PLUS;        b->signbit = PLUS;initialize_bignum(&row);initialize_bignum(&tmp);c->lastdigit = a->lastdigit;for (i=a->lastdigit; i>=0; i--) {digit_shift(&row,1);row.digits[0] = a->digits[i];c->digits[i] = 0;while (compare_bignum(&row,b) != PLUS) {c->digits[i] ++;subtract_bignum(&row,b,&tmp);row = tmp;}}zero_justify(c);a->signbit = asign;b->signbit = bsign;}bignum zero;int main(){//initint_to_bignum(0,&zero);int a[g_order * g_order] = {41, 18467, 6334, 26500,   19169, 15724, 11478, 29358,   26962, 24464, 5705, 28145,   23281, 16827, 9961, 491};//double a[g_order * g_order] = {1,2,3,4,3,2,1,3,2};fractor f[g_order * g_order];for (int i=0; i<g_order * g_order; i++) {int_to_bignum(a[i], &f[i].top);int_to_bignum(1, &f[i].down);}display(f);bignum r = det(f, 0);printf("Result = ", r);print_bignum(&r);getchar();} void exchange(fractor* matrix, int startRowColumn, int index) {for (int i=startRowColumn; i<g_order; i++) {bignum tt;copy_bignum(&tt, &M(matrix, i, index).top);bignum td;copy_bignum(&td, &M(matrix, i, index).down);M(matrix, i, index).top = M(matrix, i, startRowColumn).top;M(matrix, i, index).down = M(matrix, i, startRowColumn).down;M(matrix, i, startRowColumn).top = tt;M(matrix, i, startRowColumn).down = td;}}bignum det(fractor* matrix, int startRowColumn) {fractor result;int_to_bignum(1, &result.top);int_to_bignum(1, &result.down);//处理第i行for (int i=startRowColumn; i<g_order; i++) {if (i == g_order - 1) {bignum top, down;multiply_bignum(&result.top, &M(matrix, g_order - 1, g_order - 1).top, &top);multiply_bignum(&result.down, &M(matrix, g_order - 1, g_order - 1).down, &down);copy_bignum(&result.top, &top);copy_bignum(&result.down, &down);break;}int index = findNoneZero(matrix, i);if (index == -1) {return zero;}if (index != i) {printf("Exchanged!\n");exchange(matrix, startRowColumn, index);//交换了两列，变号allSign *= -1;}//现在M(i, i) != 0了// display(matrix);bignum top, down;multiply_bignum(&result.top, &M(matrix, i, i).top, &top);multiply_bignum(&result.down, &M(matrix, i, i).down, &down);copy_bignum(&result.top, &top);copy_bignum(&result.down, &down);// printf("Result = %lf %lf\n", result.top, result.down);fractor a = M(matrix, i, i);bignum t, d;copy_bignum(&t, &a.top);copy_bignum(&d, &a.down);// printf("a = %lf/%lf\n", t, d);for (int j = i + 1; j < g_order; j++) {bignum b, c;copy_bignum(&b, &M(matrix, i, j).top);copy_bignum(&c, &M(matrix, i, j).down);// printf("%lf %lf \n", b, c);bignum top, down;multiply_bignum(&b, &d, &top);multiply_bignum(&c, &t, &down);copy_bignum(&M(matrix, i, j).top, &top);copy_bignum(&M(matrix, i, j).down, &down);//printf("after = %lf/%lf\n", M(matrix, i, j).top, M(matrix, i, j).down);} display(matrix);// printf("extract a = %f\n", (double)a.top / (double) a.down);//然后把矩阵内其它元素变换for (int j = i + 1; j<g_order; j++) {for (int k = i + 1; k<g_order; k++) {//M(matrix, j, k) = M(matrix, j, k) //         - M(matrix, j, i) * M(matrix, i, k);bignum a, b;copy_bignum(&a, &M(matrix, j, k).top);copy_bignum(&b, &M(matrix, j, k).down);bignum c,d;multiply_bignum(&M(matrix, j, i).top, &M(matrix, i, k).top, &c);multiply_bignum(&M(matrix, j, i).down, &M(matrix, i, k).down, &d);bignum rt, rt1, rt2;multiply_bignum(&a, &d, &rt1);multiply_bignum(&b, &c, &rt2);subtract_bignum(&rt1, &rt2, &rt);bignum rd;multiply_bignum(&b, &d, &rd);copy_bignum(&M(matrix, j, k).top, &rt);copy_bignum(&M(matrix, j, k).down, &rd);//printf(" set (%d, %d) = %f\n", j, k, M(matrix, j ,k));}}display(matrix);printf("\n");}//printf("Result = %lf %lf\n", result.top, result.down);bignum r;divide_bignum(&result.top, &result.down, &r);if (allSign == -1) {r.signbit = r.signbit == PLUS ? MINUS : PLUS;}return r;}int findNoneZero(fractor* matrix, int row) {for (int i=row; i<g_order; i++) {if (compare_bignum(&M(matrix, row, i).top, &zero) != 0) {return i;}}return -1;}