一元三次方程求解——新的尝试：可变参数

/*
* 编写人    黄铎彦
* 编写日期  2023/1/11
* 程序功能  求解福州大学程序设计练习与竞赛系统 lq1086 题，一元三次方程求解
*           蓝桥杯原题链接 http://lx.lanqiao.cn/problem.page?gpid=T85
* 修改日期  2023/1/16
* 修改原因  函数 dichotomy 中有多处相同的函数调用表达式，
            将它们用变量存起来，减少调用次数
*/
/*
* 测试输入
-2 7 -5 -1

* 测试输出
-0.16 1.32 2.34
*/
#define     _CRT_SECURE_NO_WARNINGS
#define     MIN     -100    // 根最小的可能值
#define     MAX     100     // 根最大的可能值
#define     EPS     0.01    // 根的精度
#include    <stdio.h>
#include    <math.h>
#include    <stdarg.h>

int     solve_quadratic(const double a, const double b, const double c,
    double* const x1, double* const x2);
int     solve_cubic(const double a, const double b, const double c, 
    const double d, double* const x1, double* const x2, double* const x3);
double  cubic_func(const double x, ...);
int     dichotomy(double left,
    double right, double(*func)(double), double* const solution);
int     cubic_have3real(const double a, 
    const double b, const double c, const double d);

int main(void)
{
    double a; // 三次项系数
    double b; // 二次项系数
    double c; // 一次项系数
    double d; // 常数项
    scanf("%lf%lf%lf%lf", &a, &b, &c, &d);
    double x1, x2, x3;
    int err = solve_cubic(a, b, c, d, &x1, &x2, &x3);
    if (err)    puts("求解过程中出现了错误。");
    else        printf("%.2f %.2f %.2f", x1, x2, x3);
    return 0;
}

/// <summary>
/// 求一元二次方程的根
/// <para>依赖：math.h、stdio.h</para>
/// <para>保证第一个根不大于第二个根</para>
/// </summary>
/// <param name="a">二次项系数</param>
/// <param name="b">一次项系数</param>
/// <param name="c">常数项</param>
/// <param name="x1">参数输出：第一个根</param>
/// <param name="x2">参数输出：第二个根</param>
/// <returns>如果有实根，返回 0；否则返回 1。</returns>
int solve_quadratic(const double a, const double b, const double c,
    double* const x1, double* const x2)
{
    double delta = fma(b, b, -4 * a * c);
    if (delta < 0)
    {
        puts("判别式为负，没有实根！");
        return 1;
    }

    double mid = b / a / -2;
    if (delta == 0) *x1 = *x2 = mid;
    else
    {
        double offset = fabs(sqrt(delta) / a / 2); // 必须加绝对值！
        *x1 = mid - offset;
        *x2 = mid + offset;
    }
    return 0;
}

/// <summary>
/// 求一元三次方程的根。注意：参数必须保证有三个不同的实根。
/// </summary>
/// <para>依赖：头文件 stdio.h；函数 solve_quadratic、
/// cubic_func、cubic_have3real 和 dichotomy；宏 MAX 和 MIN</para>
/// <para>保证第一个根小于第二个根，且第二个根小于第三个根</para>
/// <param name="a">三次项系数</param>
/// <param name="b">二次项系数</param>
/// <param name="c">一次项系数</param>
/// <param name="d">常数项</param>
/// <param name="x1">参数输出：第一个根</param>
/// <param name="x2">参数输出：第二个根</param>
/// <param name="x3">参数输出：第三个根</param>
/// <returns>
/// 如果方程没有三个相异的实根，或者实根超出了最小或最大的可能值，
/// 返回 1；否则返回 0。
/// </returns>
int solve_cubic(const double a, const double b, const double c, 
    const double d, double* const x1, double* const x2, double* const x3)
{
    /* 计算判别式来检查是否有三个相异的实根 */
    if (!cubic_have3real(a, b, c, d))
    {
        puts("根据判别式，不存在三个相异的实根！");
        return 1;
    }

    /* 根据一阶导数求出极值点 */
    double m, n; // 极值点，其中 m < n
    int err = solve_quadratic(3 * a, 2 * b, c, &m, &n);
    if (err)
    {
        puts("求极值点时出错！");
        return 1;
    }
    if (m < MIN || n > MAX)
    {
        puts("极值超出了根最小或最大的可能值！");
        return 1;
    }

    /* 根据二分法，在左端点、两个极值点和右端点形成的三个区间中求根 */
    cubic_func(0, a, b, c, d); // 初始化三次函数
    err =   dichotomy(MIN, m, cubic_func, x1)
        ||  dichotomy(m, n, cubic_func, x2)
        ||  dichotomy(n, MAX, cubic_func, x3);
    if (err)
    {
        puts("二分法求根失败！");
        return 1;
    }
    return 0;
}

/// <summary>
/// 求三次函数的函数值
/// <para>本函数的用法是：</para>
/// <para>首次使用本函数，需要对自变量的三次项、二次项和一次项系数，
/// 以及常数项进行初始化，方法是，将它们按照顺序，作为第 2 个至第 5 个参数，
/// 传入函数内。这样，以后调用本函数，只需要一个参数作自变量就可以了。</para>
/// <para>注意：初始化时传递的额外参数必须是 double 型！</para>
/// <para>依赖：stdarg.h 和 math.h</para>
/// </summary>
/// <param name="x">自变量</param>
/// <param name="...">
/// 初始化参数（double 型）：三次项、二次项和一次项系数，以及常数项
/// </param>
/// <returns>关于自变量的三次函数值</returns>
double cubic_func(const double x, ...)
{
    static int first_call = 1;
    static double a; // 三次项系数
    static double b; // 二次项系数
    static double c; // 一次项系数
    static double d; // 常数项
    if (first_call)
    {
        first_call = 0;
        va_list ap;
        va_start(ap, x);
        a = va_arg(ap, double);
        b = va_arg(ap, double);
        c = va_arg(ap, double);
        d = va_arg(ap, double);
        va_end(ap);
    }
    return fma(x, fma(x, fma(x, a, b), c), d);
}

/// <summary>
/// 利用二分法求零点
/// <para>依赖：头文件 stdio.h 和宏 EPS</para>
/// </summary>
/// <param name="left">左端点</param>
/// <param name="right">右端点</param>
/// <param name="func">函数</param>
/// <param name="solution">参数输出：解</param>
/// <returns>
/// 如果区间两端点的函数值同号，或者左端点比右端点大，
/// 则无法求解，返回 1；否则返回 0。
/// </returns>
int dichotomy(double left, 
    double right, double(*func)(double), double* const solution)
{
    double f_left = func(left), f_right = func(right);
    if (left >= right || f_left * f_right > 0)
    {
        puts("端点非法，无法求零点！");
        return 1;
    }

    if      (f_left     == 0)   *solution = left;
    else if (f_right    == 0)   *solution = right;
    else
    {
        double mid = (left + right) / 2, f_mid = func(mid);
        for (; right - left >= EPS && f_mid != 0; 
            mid = (left + right) / 2, f_mid = func(mid))
        {
            if      (f_left * f_mid     < 0)    right   = mid;
            else if (f_mid  * f_right   < 0)    left    = mid;
        }
        *solution = mid;
    }
    return 0;
}

/// <summary>
/// 通过判别式，判断由函数参数确定的一元三次方程是否有三个不同的实根
/// <para>依赖：math.h</para>
/// </summary>
/// <param name="a">三次项系数</param>
/// <param name="b">二次项系数</param>
/// <param name="c">一次项系数</param>
/// <param name="d">常数项</param>
/// <returns>如果有三个不同的实根，返回 1；否则返回 0。</returns>
int cubic_have3real(const double a, 
    const double b, const double c, const double d)
{
    double p = c / a - pow(b, 2) / pow(a, 2) / 3;
    double q = d / a - b * c / pow(a, 2) / 3 + 2 * pow(b, 3) / pow(a, 3) / 27;
    return pow(q, 2) / 4 + pow(p, 3) / 27 < 0;
}