单精度环境下四次方程正实根求解算法故障排查与适配

我在微控制器上使用从Stack Overflow获取的单精度三次、四次方程求解算法（原双精度版本可正常运行），但单精度环境下，当cubic函数中q=0时会返回0，导致无法得到正确的实根。比如以下测试案例中，Matlab能正确求出正实根31.7814，但当前代码无法得到该结果。

原算法代码 #

float32_t cubic(float32_t b,float32_t c,float32_t d)
{
    float32_t p=c-b*b/3.0f;
    float32_t q=2.0f*b*b*b/27.0f-b*c/3.0f+d;

    if(fabsf(p)==0.0f) return powf(q,1.0f/3.0f);
    if(fabsf(q)==0.0f){
        PRINTF(INFO, "q=0 %f", p);
        return 0.0f; // TODO
    }

    float32_t t=sqrtf(fabsf(p)/3.0f);
    float32_t g=1.5f*q/(p*t);
    if(p>0.0f)
    return -2.0f*t*sinhf(asinhf(g)/3.0f)-b/3.0f;


    if(4.0f*p*p*p+27.0f*q*q<0.0f)
    return 2.0f*t*cosf(acosf(g)/3.0f)-b/3.0f;

    if(q>0.0f)
    return -2.0f*t*coshf(acoshf(-g)/3.0f)-b/3.0f;

    return 2.0f*t*coshf(acoshf(g)/3.0f)-b/3.0f;
}

int quartic(float32_t b,float32_t c,float32_t d,float32_t e,float32_t* ans)
{

    float32_t p=c-0.375f*b*b;
    float32_t q=0.125f*b*b*b-0.5f*b*c+d;
    float32_t m=cubic(p,0.25f*p*p+0.01171875f*b*b*b*b-e+0.25f*b*d-0.0625f*b*b*c,-0.125f*q*q);
    if(fabsf(q)==0.0f)
    {
        if(m<0.0f) return 0;
        int nroots=0;
        float32_t sqrt_2m=sqrtf(2.0f*m);
        if(-m-p>0.0f)
        {
            float32_t delta=sqrtf(2.0f*(-m-p));
            ans[nroots++]=-0.25f*b+0.5f*(sqrt_2m-delta);
            ans[nroots++]=-0.25f*b-0.5f*(sqrt_2m-delta);
            ans[nroots++]=-0.25f*b+0.5f*(sqrt_2m+delta);
            ans[nroots++]=-0.25f*b-0.5f*(sqrt_2m+delta);
        }

        if(fabsf(-m-p)==0.0f)
        {
            ans[nroots++]=-0.25f*b-0.5f*sqrt_2m;
            ans[nroots++]=-0.25f*b+0.5f*sqrt_2m;
        }

        return nroots;
    }

    if(m<0.0f) return 0;
    float32_t sqrt_2m=sqrtf(2.0f*m);
    int nroots=0;
    if(fabsf(-m-p+q/sqrt_2m)>=0.0f)
    {
        float32_t delta=sqrtf(2.0f*(-m-p+q/sqrt_2m));
        ans[nroots++]=0.5f*(-sqrt_2m+delta)-0.25f*b;
        ans[nroots++]=0.5f*(-sqrt_2m-delta)-0.25f*b;
    }

    if(fabsf(-m-p-q/sqrt_2m)>=0.0f)
    {
        float32_t delta=sqrtf(2.0f*(-m-p-q/sqrt_2m));
        ans[nroots++]=0.5f*(sqrt_2m+delta)-0.25f*b;
        ans[nroots++]=0.5f*(sqrt_2m-delta)-0.25f*b;
    }

    return nroots;
}

测试案例（Matlab代码） #

C1=53.3456154
C2=1729.59448
C3=54973.8164
    
C4=56.3456192
C5=1729.5946
C6=54973.8242
        
ans=single(roots([C5 (-C1+2*C4-C6) (3*C2-3*C5) (C1-2*C3+C6) -C2]))
r=r(r==conj(r));
r=r(r>0)

Matlab运行结果 #

ans =

   31.7814 +       0i
   -0.0000 +  1.0001i
   -0.0000 -  1.0001i
   -0.0315 +       0i

r = 31.7814

修复方案 #

问题出在cubic函数的q=0分支，原代码直接返回0是错误的。当q=0时，三次方程退化为：
$$(x + \frac{b}{3})^3 + p(x + \frac{b}{3}) = 0$$
提取公因式后可得实根为$x = -\frac{b}{3}$（当$p≥0$时为唯一实根；当$p<0$时，额外有两个实根$x = -\frac{b}{3} ± \sqrt{-p}$）。

修改cubic函数中的q=0分支逻辑，同时考虑单精度浮点误差，用小epsilon替代直接判断等于0：

if(fabsf(q) < 1e-6f){ // 用epsilon避免浮点误差误判
    float32_t x0 = -b / 3.0f;
    // 优先返回正实根，匹配Matlab的需求
    if(x0 > 0.0f) {
        return x0;
    }
    // 当p<0时，存在另外两个实根
    if(p < -1e-6f) {
        float32_t sqrt_neg_p = sqrtf(-p);
        float32_t x1 = x0 + sqrt_neg_p;
        float32_t x2 = x0 - sqrt_neg_p;
        
        // 返回最大的正根，若存在多个正根
        if(x1 > 0.0f && x2 > 0.0f) {
            return fmaxf(x1, x2);
        } else if(x1 > 0.0f) {
            return x1;
        } else if(x2 > 0.0f) {
            return x2;
        }
    }
    // 无正根时返回x0（可能非正）
    return x0;
}

修改后，cubic函数能正确返回q=0场景下的实根，进而quartic函数可以得到测试案例中的正实根31.7814。

内容的提问来源于stack exchange，提问作者luistope