计算行列式的值

最新推荐文章于 2024-03-19 21:23:43 发布
最新推荐文章于 2024-03-19 21:23:43 发布
阅读量3.1k 收藏 2
点赞数 1
分类专栏： 理论 文章标签： C 行列式 计算 矩阵
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/lff0305/article/details/51344355
版权

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

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

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

基本原理就是变成上三角为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()
{
	//init
	double 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"


#define	MAXDIGITS	1024		/* maximum length bignum */ 

#define PLUS		1		/* 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 4
struct 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()
{
	//init
	
	int_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;
}


lff0305
关注
专栏目录
利用Numpy计算行列式
howard2005的专栏
05-21 3157
行列式的定义、行列式的性质、行列式的性质、用数学方法计算行列式、利用Numpy工具计算行列式
求行列式的值求行列式的值求行列式的值
11-23
求行列式的值 求行列式的值 求行列式的值 求行列式的值
行列式求值
qq_41730491的博客
11-17 879
行列式(逆序数法)求值
求行列式的值
qq_39123241的博客
04-30 3982
C语言求行列式的值如何用C语言求行列式子的值呢值呢？废话少说，上代码首先定义一个函数，求行列式去掉第一行某个元素的余子式然后定义一个递归函数，求一个行列式的值然后是我们的main函数测试 如何用C语言求行列式子的值呢值呢？ 肯定要用递归，肯定要用数组 废话少说，上代码 首先定义一个函数，求行列式去掉第一行某个元素的余子式 用一个二维数组来表示这个矩阵，但是参数传递的时候我们不能确定这个矩阵的行数列...
用C语言编程递归求行列式的值
Purpleendurer@CSDN
07-10 3327
/************************************************************************功能: 用递归方法求行列式的值作者: 黄志斌 于 广西河池在自学《高等数学(二)》时，做行列式习题所得结果有些与教材所附答案不一致，为了验证究竟答案的正确性，编写了这个程序。其实用MS Office中Excel也可以验证，不过自己 DIY  感觉更舒服。
c语言编写程序计算行列式值,求行列式的值,用C语言怎么写啊? 如何求行列式的值...
weixin_29817227的博客
05-17 2087
导航：网站首页 >求行列式的值,用C语言怎么写啊? 如何求行列式的值求行列式的值,用C语言怎么写啊? 如何求行列式的值相关问题:匿名网友:行列式计算公式知道吧，给你个源码参考，大一时写的。#include #include #define MAX 9 /*定义最大为9阶行列式*/int Fun(int n, int a[MAX][MAX]...
C#使用加边法计算行列式的值
01-01
本文实例讲述了C#使用加边法计算行列式的值。分享给大家供大家参考。具体如下： 1.函数 行列式的值等于其第一行各元素乘以各自对应的代数余子式之积的和。 （注：本代码仅提供一种思路，并不代表最优解） /// /// ...
矩阵行列式计算_计算行列式_矩阵行列式计算_
10-03
当矩阵阶为2时，直接返回行列式值；否则，根据递归公式进行计算。 5. **输出结果**：最后，调用`determinant`函数并打印结果。 在实现过程中，还需要注意内存管理，确保正确地分配和释放内存。同时，为了提高效率...
hanglieshi.zip_行列式计算_计算行列式
09-21
行列式是线性代数中的一个重要概念，它在矩阵理论、几何变换、...用户可以通过编译和运行这些源代码，输入自定义的矩阵来获取相应的行列式值。这在没有图形用户界面的环境下特别有用，例如在命令行界面或嵌入式系统中。
C 计算行列式的值
04-06
C语言，自己编的。手动输入行列式，然后计算行列式的值。
行列式计算算法
07-05
本代码可以正确快速的计算出行列式，方便移植，可以运行
行列式的计算方法(含四种，看完就会！)
热门推荐
weixin_46664967的博客
02-05 60万+
行列式的计算 前言 一、对角线法 二、代数余子式法 三、等价转化法 四、逆序数法 总结 本文主要讲述行列式的求解方法，所以本文侧重于方法的讲解，而并非推导。主要思路为从三阶行列式举例，再过渡到高阶行列式的通用方法
c语言编程计算行列式的值,新手作品：行列式计算C语言版
weixin_36303045的博客
05-20 98
该楼层疑似违规已被系统折叠隐藏此楼查看此楼对话 ControlHeightDecrease Shift+Up Arrow 向上调整选定的控件或对话一个对话单位对话 ControlHeightIncrease Shift+Down Arrow 向下调整选定的控件或对话一个对话单位对话 ControlMoveDown Down Arr...
《线性代数：行列式》：常见行列式的值
逐梦者
08-29 4726
1、主对角行列式 2、副对角行列式 3、拉普拉斯展开式 设 A 是 m 阶矩阵， B 是 n 阶矩阵，则 主对角线 副对角线 4、范德蒙行列式 ...
求矩阵行列式的值（方阵）——C
berrr_的博客
05-17 5166
#include <stdio.h> #include <stdlib.h> #include <math.h> #define MAX 100 //宏定义，MAX的值为100,而且必须大写 double D(double arr[][MAX],int n); //定义一个函数来求子矩阵行列式的值； //二位数组定义必须要有列的初始化，行可以空着 void Display(double arr[][MAX],int n); //输出矩阵的函数 int main() .
线性代数|行列式定义及其值
最新发布
分享测绘及相关领域的知识
03-19 914
设有n2n^2n2个数，排成nnn行nnn列的数表a11a12⋯a1na21a22⋯a2n⋮⋮⋮an1an2⋯anna11​a21​⋮an1​​a12​a22​⋮an2​​⋯⋯⋯​a1n​a2n​⋮ann​​作出表中位于不同行不同列的nnn个数的乘积，并冠以符号−1t(-1)^t−1t，得到形如−1ta1p1a2p。
C语言实现计算行列式值的递归程序
"该代码是用于计算行列式的值的一个C语言程序。程序通过递归方式实现，用户可以输入阶数及行列式的元素，程序会输出行列式的值。" 在这个程序中，主要涉及到以下几个知识点： 1. **行列式**: 行列式是线性代数中的...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值