0

I am using these algorithms on a microcontroller:

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

From: Specialised algorithm to find positive real solutions to quartic equations? (I can't comment on this thread because I haven't enough Stackoverflow privileges)

Is the same but with single precision instead of doubles. With doubles it seems that works. With singles, i have some points in which q=0 in the cubic function, so it returns m=0 for the quartic function and the result is not a real root. I get zeros.

I need to implement Matlab root method for real positive roots.

This point gives problems with the above code but has solution with roots in 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)

In Matlab the result is ok

ans =

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

r = 31.7814
Jason Aller
  • 3,541
  • 28
  • 38
  • 38
luistope
  • 1
  • 2
  • 2
    `fabsf(-m-p)==0.0f` Do you really mean to compare w zero? This is typically not what you want to do w floating point numbers. Typically you have a tolerance for the relative error based on machine epsilon. – wcochran Nov 23 '22 at 17:35
  • Please provide an example call of the `quartic` function, the actual result, and the expected result. – Ian Abbott Nov 23 '22 at 17:46
  • Are you committed to this algorithm? If you have a good eigen solver routine you can build the companion matrix C for the polynomial and find the eigenvalues of C. http://web.mit.edu/18.06/www/Spring17/Eigenvalue-Polynomials.pdf – wcochran Nov 23 '22 at 17:57
  • @Ian Abbott `b = (-31.7516479); p = (-378.062683); q = (-4033.1228); c = p + 0.375*b*b; d = q - 0.125*b^3+0.5*b*c; e = (-1729.59448/1729.5946); %e = -1; % Q in cubic q_c = (2/27)*(p)^3 - (p)*(0.25*p^2+0.01171875*b^4-e+0.25*b*d-0.0625*c*b^2)/3 + (-0.125*q^2); ` It returns: q_c = -6.9952e-04 If I use single for all the variables in Matlab q_c = 0. It is the same problem inside the microcontroller. Comparison == 0.0f gives true and returns 0.0f. Not correct. With doubles, I do not have this problem (I could have the same problem with lower values). – luistope Nov 24 '22 at 11:17
  • @wcochran I avoid iterative methods and prefer a close-form method to have a deterministic calculation time. How to avoid comparison with zero? It's the same to compare < than an epsilon. If the result is near 0 I get the same result. Or there is another solution. I need to make comparisons to distinguish between different cases inside the solutions (Three real roots, one real,...). – luistope Nov 24 '22 at 11:23

0 Answers0