Comparing a double against zero

Question

I'm new to Java and I've been trying to implement an algorithm for finding the roots of a cubical equation. The problem arises when I calculate the discriminant and try to check where it falls relative to zero.

If you run it and enter the numbers "1 -5 8 -4", the output is as follows:

1 -5 8 -4
p=-0.333333, q=0.074074
disc1=0.001372, disc2=-0.001372
discriminant=0.00000000000000001236
Discriminant is greater than zero.

I know the problem arises because the calculations with doubles are not precise. Normally the discriminant should be 0, but it ends up being something like 0.00000000000000001236.

My question is, what is the best way to avoid this? Should I check if the number falls between an epsilon neighborhood of zero? Or is there a better and more precise way?

Thank you in advance for your answers.

import java.util.Scanner;

class Cubical {
    public static void main(String[] args) {
        // Declare the variables.
        double a, b, c, d, p, q, gamma, discriminant;

        Scanner userInput = new Scanner(System.in);
        a = userInput.nextDouble();
        b = userInput.nextDouble();
        c = userInput.nextDouble();     
        d = userInput.nextDouble();

        // Calculate p and q.
        p = (3*a*c - b*b) / (3*a*a);
        q = (2*b*b*b) / (27*a*a*a) - (b*c) / (3*a*a) + d/a;

        // Calculate the discriminant.
        discriminant = (q/2)*(q/2) + (p/3)*(p/3)*(p/3);

        // Just to see the values.
        System.out.printf("p=%f, q=%f\ndisc1=%f, disc2=%f\ndiscriminant=%.20f\n", p, q, (q/2)*(q/2), (p/3)*(p/3)*(p/3), (q/2)*(q/2) + (p/3)*(p/3)*(p/3));

        if (discriminant > 0) {
            System.out.println("Discriminant is greater than zero.");
        }
        if (discriminant == 0) {
            System.out.println("Discriminant is equal to zero.");
        }
        if (discriminant < 0) {
            System.out.println("Discriminant is less than zero.");
        }
    }
}

Answer 1

The simplest epsilon check is

if(Math.abs(value) < ERROR)

a more complex one is proportional to the value

if(Math.abs(value) < ERROR_FACTOR * Math.max(Math.abs(a), Math.abs(b)))

---

In your specific case you can:

if (discriminant > ERROR) {
    System.out.println("Discriminant is greater than zero.");
} else if (discriminant < -ERROR) {
    System.out.println("Discriminant is less than zero.");
} else {
    System.out.println("Discriminant is equal to zero.");
}

Answer 2

Here's solution that is precise when the input values are integers, though it is probably not the most practical.

It will probably also work fine on input values that have a finite binary representation (eg. 0.125 does, but 0.1 doesn't).

The trick: Remove all divisions from the intermediate results and only divide once at the end. This is done by keeping track of all the (partial) numerators and denominators. If the discriminant should be 0 then it's numerator will be 0. No round-off error here as long as values at intermediate additions are within a magnitude of ~2^45 from each other (which is usually the case).

// Calculate p and q.
double pn = 3 * a * c - b * b;
double pd = 3 * a * a;

double qn1 = 2 * b * b * b;
double qd1 = 27 * a * a * a;
double qn2 = b * c;
double qn3 = qn1 * pd - qn2 * qd1;
double qd3 = qd1 * pd;
double qn = qn3 * a + d * qd3;
double qd = qd3 * a;

// Calculate the discriminant.
double dn1 = qn * qn;
double dd1 = 4 * qd * qd;
double dn2 = pn * pn * pn;
double dd2 = 27 * pd * pd * pd;

double dn = dn1 * dd2 + dn2 * dd1;
double dd = dd1 * dd2;
discriminant = dn / dd;

(only checked on the provided input values, so tell me if something's wrong)

Answer 3

Should I check if the number falls between an epsilon neighborhood of zero?

Exactly

Answer 4

maybe BigDecimal is worth a look at...

http://download.oracle.com/javase/1.4.2/docs/api/java/math/BigDecimal.html

you can secify the round mode in the divide-operation