一元三次/一元二次方程求数值解C++编程实践

一元三次方程

        一元三次方程求解目前最简明的、最直观的方法是：盛金公式法。这是80年代中国中学数学教师范盛金的研究成果，体现了体现了数学的有序、对称、和谐与简洁美。

        盛金公式法详细展开请参见一元三次方程求解(求根) - 盛金公式法-CSDN博客文章浏览阅读2.8w次，点赞26次，收藏124次。一元三次方程求解(求根) - 盛金公式法一、引言 只含有一个未知数（即“元”），并且未知数的最高次数为3（即“次”）的整式方程叫做一元三次方程（英文名：cubic equation in one unknown）。一元三次方程的标准形式（即所有一元三次方程经整理都能得到的形式）是ax3+bx2+cx+d=0（a，b，c，d为常数，x..._盛金公式https://blog.csdn.net/u012912039/article/details/101363323

C++编程实践

        一元三次方程求数值解涉及判别式Δ与0比大小，求三角反三角以、求平方根和求立方根，因此C++编程涉及两个细节：

        1、考虑浮点运算精度问题，Δ与0比大小需要设一个阈值eps啊（用Δ>eps取代Δ>0，用Δ>-eps并且Δ<eps取代Δ=0，用Δ<-eps取代Δ<0；

        2、负数求立方根的处理，先对负数的相反数求立方根再取相反数，即-pow(-x, 1./3.)

C++代码模板

#include <cmath>

#define eps 1e-12

double cubic_root(double x) {
    return x<0. ? -pow(-x, 1./3.) : pow(x, 1./3.);
}

int solve_cubic_equation(double a, double b, double c, double d, double (&x)[3]) {
    // 一元三次方程求解，返回不同实根的数量，a!=0
    double A = b*b - 3*a*c, B = b*c - 9*a*d, C = c*c - 3*b*d, delta = B*B - 4.*A*C;
    // https://blog.csdn.net/u012912039/article/details/101363323
    // A=B=0时，方程有三重实根-c/b
    // delta > 0 时，方程有一个实根和一对共轭虚根
    // delta = 0 时，方程有有两个不同实根，其中一个是两重根
    // delta < 0 时，方程有三个不同实根
    if (A == 0. && B == 0.) {
        x[0] = -c/b;
        return 1;
    } else if (delta > 0.) {
        delta = sqrt(delta);
        double y1 = A*b + 1.5*a*(delta-B), y2 = A*b - 1.5*a*(delta+B);
        x[0] = -(cubic_root(y1) + cubic_root(y2) + b)/3./a;
        return 1;
    } else if (delta > -eps) {
        x[0] = B/A - b/a; x[1] = -B/2./A;
        return 2;
    } 
    double sA = sqrt(A), t = acos((A*b-1.5*a*B)/A/sA)/3., s = sqrt(3)*sin(t); t = cos(t); a *= -3.;
    x[0] = (b + 2*sA*t) / a; x[1] =  (b - sA*(t + s)) / a; x[2] = (b - sA*(t - s)) / a;
    return 3;
}

double find_a_root(double a, double b, double c, double d) {
    // 求出一元三次方程的一个实根，a!=0
    double A = b*b - 3*a*c, B = b*c - 9*a*d, C = c*c - 3*b*d, delta = B*B - 4.*A*C;
    if (A == 0. && B == 0.) return -c/b;
    if (delta > 0.) {
        delta = sqrt(delta);
        double y1 = A*b + 1.5*a*(delta-B), y2 = A*b - 1.5*a*(delta+B);
        return -(cubic_root(y1) + cubic_root(y2) + b)/3./a;
    } else if (delta > -eps) return -B/2./A;
    double sA = sqrt(A), t = acos((A*b-1.5*a*B)/A/sA)/3.;
    return -(b + 2*sA*cos(t)) / 3. / a;
}

int solve_quadratic_equation(double a, double b, double c, double *x) {
    // 解一元二次方程（a可以为0）并返回实根数量
    if (a == 0.) {
        if (b == 0.) return 0;
        *x = -c/b;
        return 1;
    }
    double delta = b*b - 4.*a*c;
    if (delta < 0.) {
        if (delta < -eps) return 0;
        delta = 0.;
    }
    a *= 2; delta = sqrt(delta);
    if (delta == 0. || delta < a*eps) {
        *x = -b/a;
        return 1;
    }
    *x = (-b-delta)/a; *(++x) = (delta-b)/a;
    return 2;
}