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.

As a side note you can definitely look into caching/memoization of factorial results to speed up the program, because that is being computed repeatedly. But that is not a bug, since you are getting a result. None of these calculations should overflow `double` at n=14. It would be helpful if you can share more about the expected vs actual results for each n value. — shriakhilc, Dec 27 '21 at 03:21
Put any necessary information in your question. And I doubt that the Wolfram article tells us what numbers you are actually seeing. — tgdavies, Dec 27 '21 at 03:57
Sorry, 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. — Babai08, Dec 27 '21 at 04:10
I think that large intermediate values are pushing precision you need in the final result out of the mantissa. — tgdavies, Dec 27 '21 at 05:36
So do you think its the size of the factorial because yeah 14! is big — Babai08, Dec 27 '21 at 05:38

tgdavies · Accepted Answer · 2021-12-28T05:14:33.477

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.

Thank you so much. Although it still diverges past there, I'm getting 15: -23.639285714285716 16: -2245.162672376643 17: -14389.589377452612 18: -147997.97879116068 19: 7.928308554807248E8 20: 2.275498450290969E10 and if I go any further it halts. — Babai08, Dec 27 '21 at 06:14
v^n was losing accuracy, I've changed that to use bigdecimal — tgdavies, Dec 28 '21 at 00:00

Issues with accuracy with Bernoulli Number generator in java

1 Answers1