QE solve and a check if the roots satisfy the equation : using C language

Question

I am writing a code in C. Checking the roots if they satisfy the quadratic equation. In printf, while %lf outputs 0, %e returns a very small value (between 10^-15 and 10^-19). I tried to declare the variable as float or double, nothing is working when I am printing using %e. What is wrong?

Also any suggestions to improve the logic would be appreciated. I am a beginner in C. thanks.

----------------------code-------------------------------

#include <stdio.h>
#include <math.h>

struct param {
    double a[11];
} x, y, z;

double slvequad(double a1, double b1, double c1) {
    double quad_sol1, quad_sol2;
    double dcrm = b1 * b1 - 4.0 * a1 * c1;
    printf("\n ( %+0.2e, %+0.2e, %+0.2e)", a1, b1, c1);

    if (a1 != 0 && dcrm >= 0) {
        quad_sol1 = (-b1 + sqrt(dcrm)) / (2.0 * a1);
        quad_sol2 = (-b1 - sqrt(dcrm)) / (2.0 * a1);
        printf("( %+0.2e , %+0.2e ) ", quad_sol1, quad_sol2);
        quadcheck(a1, b1, c1, quad_sol1, quad_sol2);
    } else
    if (a1 == 0 && b1 != 0) {
        quad_sol1 = quad_sol2 = -c1 / b1;
        printf("( %+0.2e  ) ", quad_sol1);
        quadcheck(a1, b1, c1, quad_sol1, quad_sol2);
    } else {
        printf(" No solution  ");
    }
}

void quadcheck(double x1, double y1, double z1, double sol1, double sol2) {

    double v1 = (x1 * (sol1) * (sol1) + y1 * sol1 + z1);
    double v2 = (x1 * (sol2) * (sol2) + y1 * sol2 + z1);

    printf("( %+0.2e , %+0.2e ) \n ", v1, v2);
}

int main() {
    struct param x = { 1e-10, 0,  1, 0, 0, 1e-35, 1,  3, 4,   1,  1.0  };
    struct param y = { 2.0,   1,  0, 1, 0, 0,     4,  5, -20, 6,  -1e-1 };
    struct param z = { 1e-10, 5, -4, 0, 0, -1e35, 1, -7, 25,  34, 0.0025 };
    int i;

    int len = *(&x.a + 1) - x.a;

    for (i = 0; i < len; i++) {
        slvequad(x.a[i], y.a[i], z.a[i]);
    }

    return 0;
}

output: In format: (a1,b1,c1)(quad_sol1, quad_sol2)(v1,v1)

Technically, v1 and v2 should print 0, but it prints very small value in a few cases (using %e) as shown in the image.

Code is _working_. %lf returns an answer _rounded_ to the nearest 0.000001. %e provides more detail. What might be wrong is the unposted code that calculated the answer. — chux - Reinstate Monica, Jul 12 '22 at 12:46
Your code snippet linked does not show the actual value of `s1` and `s2` anywhere. Is the problem that "very small" is supposed to be exactly 0? — Weather Vane, Jul 12 '22 at 12:48
Often, the exact roots cannot be represented exactly and you need to check if they are close enough to 0, or check for a sign change in the function value close to the candidate root. — Ian Abbott, Jul 12 '22 at 12:51
@WeatherVane yes should be zero. the quadcheck was to check the roots putting back in QE, hence it should be zero. — san, Jul 12 '22 at 13:22
Please see [Is floating point math broken?](http://stackoverflow.com/questions/588004/is-floating-point-math-broken) and [Why Are Floating Point Numbers Inaccurate?](https://stackoverflow.com/questions/21895756/why-are-floating-point-numbers-inaccurate) and [Is floating-point == ever OK?](https://stackoverflow.com/questions/4682889/is-floating-point-ever-ok) — Weather Vane, Jul 12 '22 at 13:51

score 0 · Accepted Answer · answered Jul 14 '22 at 18:08

The roots for x² + 4x + 1 are irrational numbers -2 +/-sqrt(3). They cannot be represented exactly in types float or double, only approximations are stored, precise to about 17 significant digits for type double. Computing the polynomial with these approximations produce very small, yet non zero numbers.

You may want to read these questions:

QE solve and a check if the roots satisfy the equation : using C language

1 Answers1