一、问题描述

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

二、CCC代码

#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);
}

三、参考文献/资料

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

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

一元二次方程高精度实数根(C语言)
文章目录一.问题描述二.推导步骤三. C C C代码四.总结五.参考文献/资料一.问题描述已知一元二次方程 a x 2 + b x + c = 0 ( a ≠ 0 ) ax^2 + ...
一元三次方程的实数根
一元三次方程的实数根一元三次方程的一般形式一元三次方程的一般形式为 ax3+bx2+cx+d=0,a≠0. ax^3+bx^2+cx+d=0,\quad a\neq 0. 令 x=y−b3ax=y ...
C语言——二分法求一元三次方程的方程根
C语言--二分法求一元三次方程的方程根 #include <stdio.h>#include <math.h>int main() {double m, n;double i ...
java求一元二次方程“ax2＋bx＋c＝0，（a≠0）”的根。要求：通过命令行参数或者Scanner类输入a、b和c的值，借助Math类计算并输出该一元二次方程的实数根。
编程题:求一元二次方程"ax2+bx+c＝0,(a≠0)"的根.要求:通过命令行参数或者Scanner类输入a.b和c的值,借助Math类计算并输出该一元二次方程的实数根. 提示: ...
C语言例题——求一元二次方程的实数根
#include "stdio.h" #include "math.h"void main() {//方程的系数int a, b, c;//△和方程的解doub ...
1205: 求一元二次方程的实数根（C）
一.题目 acm.wust.edu.cn/problem.php?id=1205&soj=0 二.分析一元二次方程有三个系数a.b.c,两个根x1.x2,以及d(德尔塔): a.b.c均为实 ...
C语言一元三次方程二分法切线法
想起了自己之前刚学习C语言的时候总喜欢做一些有意思的事,最让自己满意的就是解一元三次方程.这其中也有别的原因的,高中时候也研究过一元三次方程,但是当时绞尽脑汁也没有想出求解的办法,只是找出了三个根之间 ...
MATLAB实现一元三次方程求解/盛金公式
MATLAB实现一元三次方程求解/盛金公式一元三次方程求解中,1945年卡尔丹诺把冯塔纳的三次方程求根公式发表出来,但该公式形式比较复杂,直观性也较差.1989年范盛金对一元三次方程求解进行了深入的 ...
python中一元二次方程的判别式_一元二次方程难？根的判别式常考这四题型，若只知其一难考高分...
学习有方法,解题有技巧.对于九年级学生来说,现在的学习可以说是进入"白热化"阶段,如何在紧张的学习中脱颖而出.在同样的努力前提下,对重要考点进行全面总结,不留知识死角也是要提上日程 ...

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

文章目录

一、问题描述

二、CCC代码

三、参考文献/资料

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

最新文章

热门文章