0

I am trying to calculate the roots of a quadratic equation using the Citardauq Formula, which is a more numerically stable way to calculate those roots. However, when, for example, I enter the equation x^2+200x-0.000002=0 this program does not calculate the roots precisely. Why? I don't find any error in my code and the catastrophic cancellation should not occur here.

You can find why the Citardauq formula works here (second answer).

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

int main()
{
    double a, b, c, determinant;
    double root1, root2;

    printf("Introduce coefficients a b and c:\n");
    scanf("%lf %lf %lf", &a, &b, &c);

    determinant = b * b - 4 * a * c;

    if (0 > determinant)
    {
        printf("The equation has no real solution\n");
        return 0;
    }

    if (b > 0)
    {
        root1 = (-b - sqrt(determinant)) / (2 * a);
        root2 = (c / (a * root1));
        printf("The solutions are %.16lf and %.16lf\n", root1, root2);
    }
    else if (b < 0)
    {
        root1 = (-b + sqrt(determinant)) / (2 * a);
        root2 = (c / (a * root1));
        printf("The solutions are %.16lf and %.16lf\n", root1, root2);
    }
}
codingnight
  • 231
  • 1
  • 7
  • What result did you get? What did you expect? [This old question and its answers](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) might be interesting for you to read as well. – Some programmer dude Mar 04 '18 at 15:29
  • Aside: `else if (b<0)` might be better as `else` – Weather Vane Mar 04 '18 at 15:31
  • @Someprogrammerdude I got that one root was 0.0000000100000000, getting a loss of significance, (what i wanted to avoid). Basically I just want to code [this](https://en.wikipedia.org/wiki/Loss_of_significance#A_better_algorithm) in C. – codingnight Mar 04 '18 at 16:11

1 Answers1

1

Welcome to numerical computations. There are a few issues here:

1) As pointed by some-programmer-dude there is a problem with precise representation of floating numbers Is floating point math broken?

For 0.1 in the standard binary64 format, the representation can be written exactly as 0.1000000000000000055511151231257827021181583404541015625

2) Double precision (double) gives you only 52 bits of significant, 11 bits of exponent, and 1 sign bit. Floating point numbers in C use IEEE 754 encoding.

3) sqrt precision is also limited.

In your case, the solution is as follows: enter image description here

You can see that from precision point of view it is not easy equation.

On line calculator 1 gives solutions as:

1.0000007932831068e-8  -200.00000001

Your program is better:

Introduce coefficients a b and c:                                                                                                    

1                                                                                                                                    
200                                                                                                                                  
-0.000002                                                                                                                            
The solutions are -200.0000000100000079 i 0.0000000100000000

So one of the roots is -200.000000010000. Forget about the rest of the digits. This is exactly what one can expect since double has 15 decimal digits of precision!

sg7
  • 6,108
  • 2
  • 32
  • 40
  • Thank your for your answer, I think I understand it better now. However, the second root is 0.000000010000000 and therefore I lose 7 digits of significance, right? What I wanted to achieve coding this was to avoid this loss of significance. I have another code that calculates the roots using the normal quadratic formula, and they give the same answer. Thinking about it, both programs should give the same answer, since those are the roots, but how can I avoid the loss of significance then? – codingnight Mar 04 '18 at 16:22
  • @codingnight No, you do not lose 7 digits. In `0.0000000100000000` you also have `15` "good" digits. – sg7 Mar 04 '18 at 16:27
  • Then, in [this](https://onlinegdb.com/B1fSbjY_M) code I don't lose any digits either? – codingnight Mar 04 '18 at 16:40
  • @codingnight Yes, as per double resolution, you will always have `15` "good" digits. That can be altered by numerical instability of the algorithm but it is not your case. The numbers `79` are outside of the `15` good digits range. They are just a digital "noise". – sg7 Mar 04 '18 at 16:54