一元三次方程高精度实数根(C语言)

一、问题描述

  已知一元三次方程 a x 3 + b x 2 + c x + d = 0 ( a ≠ 0 ) ax^3 + bx^2 + cx+d = 0(a\neq0) ax3+bx2+cx+d=0(a=0),求方程的实数根。笔者用C语言实现了著名的盛金公式求解一元三次方程高精度的实数根。一元三次方程实数根求解过程会调用一元二次方程高精度实数根的求解函数,参见另外一篇博文:一元二次方程高精度实数根(C语言)

二、 C C C代码

#include <math.h>
#include <float.h>
#include <stdio.h>

#define UINT32 unsigned int
#define ERR_NO_ERROR 0x00000000
#define ERR_NAN 0x00000001
#define ERR_INF 0x00000002
#define MAX(a, b) ((a) > (b)) ? (a) : (b)
#define MIN(a, b) ((a) < (b)) ? (a) : (b)

/*************************************************
Function: is_number
Description: 判断浮点数是否为nan
Input: 浮点数x
Output: 无
Return: 若浮点数x为nan返回0,否则返回1
Author: Marc Pony(marc_pony@163.com)
*************************************************/
int is_number(float x) 
{ 
	return (x == x);
}

/*************************************************
Function: is_finite_number
Description: 判断浮点数是否为inf
Input: 浮点数x
Output: 无
Return: 若浮点数x为inf返回0,否则返回1
Author: Marc Pony(marc_pony@163.com)
*************************************************/
int is_finite_number(float x) 
{ 
	return (x >= -FLT_MAX && x <= FLT_MAX);
}

/*************************************************
Function: solve_quadratic_equation
Description: 求一元二次方程(a*x^2 + b*x + c = 0)的所有实数根
Input: 方程的系数 p = {c, b, a}
Output: 方程的所有实数根x, 实数根的个数rootCount
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
*************************************************/
UINT32 solve_quadratic_equation(float p[], float x[], int *rootCount)
{
	int i;
	float a, b, c, delta, sqrtDelta;
	const float ZERO = FLT_MIN;  // min normalized positive value(1.175494351e-38F)
	const float EPS = FLT_MIN;
	UINT32 errNo = ERR_NO_ERROR;

	*rootCount = 0;
	for (i = 0; i < 3; i++)
	{
		if (!is_number(p[i]))
		{
			errNo = ERR_NAN;
			return errNo;
		}
		if (!is_finite_number(p[i]))
		{
			errNo = ERR_INF;
			return errNo;
		}
	}

	a = p[2];
	b = p[1];
	c = p[0];
	
	if (fabs(a - 0.0) < EPS)
	{
		if (fabs(b - 0.0) > EPS)
		{
			x[0] = -c / b;
			*rootCount = 1;
		}
	}
	else
	{
		b /= a;
		c /= a;
		a = 1.0;

		delta = b * b - 4.0 * a * c;
		if (delta > ZERO)
		{
			if (fabs(c - 0.0) < EPS)	//若c = 0,由于计算误差,sqrt(b*b - 4*a*c)不等于|b|
			{
				x[0] = 0.0;
				x[1] = -b / a;
			}
			else
			{
				sqrtDelta = sqrt(delta);
				if (b > 0.0)
				{
					x[0] = (-2.0 * c) / (b + sqrtDelta);	//避免两个很接近的数相减,导致精度丢失
					x[1] = (-b - sqrtDelta) / (2.0 * a);
				}
				else
				{
					x[0] = (-b + sqrtDelta) / (2.0 * a);
					x[1] = (-2.0 * c) / (b - sqrtDelta);	//避免两个很接近的数相减,导致精度丢失
				}
			}
			*rootCount = 2;
		}
		else if (fabs(delta - 0.0) < EPS)
		{
			x[0] = x[1] = -b / (2.0 * a);
			*rootCount = 2;
		}
		else
		{
			*rootCount = 0;
		}
	}
	return errNo;
}

/*************************************************
Function: solve_cubic_equation
Description: 盛金公式求一元三次方程(a*x^3 + b*x^2 + c*x + d = 0)的所有实数根
			 A = b * b - 3.0 * a * c;
			 B = b * c - 9.0 * a * d;
			 C = c * c - 3.0 * b * d;
			 (1)当A = B = 0时,方程有一个三重实根
			 (2)当Δ = B^2-4 * A * C > 0时,方程有一个实根和一对共轭虚根
			 (3)当Δ = B^2-4 * A * C = 0时,方程有三个实根,其中有一个两重根
			 (4)当Δ = B^2-4 * A * C < 0时,方程有三个不相等的实根
Input: 方程的系数 p = {d, c, b, a}
Output: 方程的所有实数根x,实数根的个数rootCount
Return: 错误号
Author: Marc Pony(marc_pony@163.com)
*************************************************/
UINT32 solve_cubic_equation(float p[], float x[], int *rootCount)
{
	int i;
	float a, b, c, d, A, B, C, delta;
	float Y1, Y2, Z1, Z2, K, parm[3], roots[2], theta, T;
	const float ZERO = FLT_MIN;  // min normalized positive value(1.175494351e-38F)
	const float EPS = FLT_MIN;
	const float CALCULATE_ERROR = 1.0e-7;
	UINT32 errNo = ERR_NO_ERROR;

	*rootCount = 0;
	for (i = 0; i < 4; i++)
	{
		if (!is_number(p[i]))
		{
			errNo = ERR_NAN;
			return errNo;
		}
		if (!is_finite_number(p[i]))
		{
			errNo = ERR_INF;
			return errNo;
		}
	}

	a = p[3];
	b = p[2];
	c = p[1];
	d = p[0];

	if (fabs(a - 0.0) < EPS)
	{
		parm[2] = b;
		parm[1] = c;
		parm[0] = d;

		errNo = solve_quadratic_equation(parm, x, rootCount);
	}
	else
	{
		b /= a;
		c /= a;
		d /= a;
		a = 1.0;

		A = b * b - 3.0 * a * c;
		B = b * c - 9.0 * a * d;
		C = c * c - 3.0 * b * d;

		delta = B * B - 4.0 * A * C;

		if (fabs(A - 0.0) < EPS && fabs(B - 0.0) < EPS)
		{
			x[0] = x[1] = x[2] = -b / (3.0 * a);
			*rootCount = 3;
			return errNo;
		}

		if (delta > ZERO)
		{
			parm[2] = 1.0;
			parm[1] = B;
			parm[0] = A * C;

			errNo = solve_quadratic_equation(parm, roots, rootCount);
			if (errNo != ERR_NO_ERROR)
			{
				return errNo;
			}
			Z1 = roots[0];
			Z2 = roots[1];

			Y1 = A * b + 3.0 * a * Z1;
			Y2 = A * b + 3.0 * a * Z2;

			if (Y1 < 0.0 && Y2 < 0.0)	//pow函数的底数必须为非负数,必须分类讨论
			{
				x[0] = (-b + pow(-Y1, 1.0 / 3.0) + pow(-Y2, 1.0 / 3.0)) / (3.0*a);
			}
			else if (Y1 < 0.0 && Y2 > 0.0)
			{
				x[0] = (-b + pow(-Y1, 1.0 / 3.0) - pow(Y2, 1.0 / 3.0)) / (3.0*a);
			}
			else if (Y1 > 0.0 && Y2 < 0.0)
			{
				x[0] = (-b - pow(Y1, 1.0 / 3.0) + pow(-Y2, 1.0 / 3.0)) / (3.0*a);
			}
			else
			{
				x[0] = (-b - pow(Y1, 1.0 / 3.0) - pow(Y2, 1.0 / 3.0)) / (3.0*a);
			}
			*rootCount = 1;
		}
		else if (fabs(delta - 0.0) < EPS)
		{
			if (fabs(A - 0.0) > EPS)
			{
				K = B / A;
				x[0] = -b / a + K;
				x[1] = x[2] = -0.5 * K;
				*rootCount = 3;
			}
		}
		else
		{
			if (A > 0.0)
			{
				T = (2.0 * A * b - 3.0 * a * B) / (2.0 * pow(A, 3.0 / 2.0));
				if (T > 1.0)	//由于计算误差,T的值可能略大于1(如1.0000001)
				{
					if (T < 1.0 + CALCULATE_ERROR)
					{
						T = 1.0;
					}
					else
					{
						return errNo;
					}
				}
				if (T < -1.0)
				{
					if (T > -1.0 - CALCULATE_ERROR)
					{
						T = -1.0;
					}
					else
					{
						return errNo;
					}
				}
				theta = acos(T);
				x[0] = (-b - 2.0 * sqrt(A) * cos(theta / 3.0)) / (3.0 * a);
				x[1] = (-b + sqrt(A) * (cos(theta / 3.0) + sqrt(3.0) * sin(theta / 3.0))) / (3.0 * a);
				x[2] = (-b + sqrt(A) * (cos(theta / 3.0) - sqrt(3.0) * sin(theta / 3.0))) / (3.0 * a);
				*rootCount = 3;
			}
		}
	}
	return errNo;
}


void main(void)
{
	float x[3], p[4];
	int rootCount;
	float a, b, c, d;
	UINT32 errNo = ERR_NO_ERROR;

	//一元三次方程测试
	//(1)(x - 1) * (x^2 + 1) = 0 (x^3 - x^2 + x - 1 = 0)
	a = 1;
	b = -1;
	c = 1;
	d = -1;

	//(2) (x - 1)^3 = 0 (x^3 - 3*x^2 + 3*x - 1 = 0)
	//a = 1;
	//b = -3;
	//c = 3;
	//d = -1;

	//(3) (x - 1)^2 * (x - 2) = 0 (x^3 - 4*x^2 + 5*x - 2 = 0)
	//a = 1;
	//b = -4;
	//c = 5;
	//d = -2;

	//(4) (x - 1) * (x - 2) * (x - 3) = 0 (x^3 - 6*x^2 + 11*x - 6 = 0)
	//a = 1;
	//b = -6;
	//c = 11;
	//d = -6;

	//(5) 0*x^3 + x^2 - 2*x + 1 = 0
	//a = 0;
	//b = 1;
	//c = -2;
	//d = 1;

	//(6) 0*x^3 + 0*x^2 - 2*x + 1 = 0
	//a = 0;
	//b = 0;
	//c = -2;
	//d = 1;

	//(7) 0*x^3 + 0*x^2 + 0*x + 1 = 0
	//a = 0;
	//b = 0;
	//c = 0;
	//d = 1;

	p[0] = d;
	p[1] = c;
	p[2] = b;
	p[3] = a;
	errNo = solve_cubic_equation(p, x, &rootCount);
}

三、参考文献/资料

三次方程新解法——盛金公式解题法

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值