Issues with accuracy with Bernoulli Number generator in java-CodePudding

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.

CodePudding user response：

Here's an implementation using BigDecimal which seems to give better results.

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.