Issues with accuracy with Bernoulli Number generator in java

Question

I've created some code that generates the Bernoulli Numbers based off of formula 33 on MathWorld. This is given at https://mathworld.wolfram.com/BernoulliNumber.html and should work for all integers n but it diverges from the expected results extremely quickly once it gets to n=14. I think the issue may be in the factorial code, although I have no idea.

It's pretty accurate up until 13, all odd numbers should be 0 besides 1 but the values past 14 give weird values. For instance 14 gives a number like 0.9 when it should give something around 7/6 and say 22 gives a very negative number in the order of 10^-4. The odd numbers give strange values like 15 gives around -11.

Here is all the related code

public static double bernoulliNumber2(int n) {
    double bernoulliN = 0;
    for (double k = 0D; k <= n; k++) {
        bernoulliN += sum2(k,n)/(k+1);
    }
    return bernoulliN;
}
public static double sum2(double k, int n) {
    double result = 0;

    for (double v = 0D; v <= k; v++) {
        result += Math.pow(-1, v) * MathUtils.nCr((int) k,(int) v) * Math.pow(v, n);
    }

    return result;    
}

public static double nCr(int n, int r) {
    return Factorial.factorial(n) / (Factorial.factorial(n - r) * Factorial.factorial(r));
}

public static double factorial(int n) {
    if (n == 0) return 1;
    else return (n * factorial(n-1));
}

Thank you in advance.

Answer 1

The problem here is that floating point arithmetic doesn't need to overflow to experience catastrophic loss of precision.

A floating point number has a mantissa and an exponent, where the value of the number is mantissa * 10^exponent (real floating point numbers use binary, I'm using decimal). The mantissa has limited precision.

When we add floating point numbers of different signs we can end up with a final result which has lost precision.

e.g. let's say the mantissa is 4 digits. If we add:

1.001 x 10^3 + 1.000 x 10^4 - 1.000 x 10^4

we expect to get 1.001 x 10^3. But 1.001 x 10^3 + 1.000 x 10^4 = 11.001 x 10^3, which is represented as 1.100 x 10^4, given that our mantissa has only 4 digits.

So when we subtract 1.000 x 10^4 we get 0.100 x 10^4, which is represented as 1.000 x 10^3 rather than 1.001 x 10^3.

Here's an implementation using BigDecimal which gives better results (and is far slower).

import java.math.BigDecimal;
import java.math.RoundingMode;

public class App {
    public static double bernoulliNumber2(int n) {
        BigDecimal bernoulliN = new BigDecimal(0);
        for (long k = 0; k <= n; k++) {
            bernoulliN = bernoulliN.add(sum2(k,n));
            //System.out.println("B:" + bernoulliN);
        }
        return bernoulliN.doubleValue();
    }
    public static BigDecimal sum2(long k, int n) {
        BigDecimal result = BigDecimal.ZERO;

        for (long v = 0; v <= k; v++) {
            BigDecimal vTon = BigDecimal.valueOf(v).pow(n);
            result = result.add(BigDecimal.valueOf(Math.pow(-1, v)).multiply(nCr(k,v)).multiply(vTon).divide(BigDecimal.valueOf(k + 1), 1000, RoundingMode.HALF_EVEN));
        }
        return result;
    }
    public static BigDecimal nCr(long n, long r) {
        return factorial(n).divide(factorial(n - r)).divide(factorial(r));
    }
    public static BigDecimal factorial(long n) {
        if (n == 0) return BigDecimal.ONE;
        else return factorial(n-1).multiply(BigDecimal.valueOf(n));
    }
    public static void main(String[] args) {
        for (int i = 0; i < 20; i++) {
            System.out.println(i + ": " + bernoulliNumber2(i));
        }
    }
}

Try changing the scale passed to the division in sum2 and watch the effect on the output.