java中伯努利数生成器的准确性问题

Question

我创建了一些代码，可以根据 MathWorld 上的公式 33 生成伯努利数。这是在https://mathworld.wolfram.com/BernoulliNumber.html 给出的，应该适用于所有整数 n，但一旦达到 n=14，它就会非常迅速地偏离预期结果。我认为问题可能出在阶乘代码中，虽然我不知道。

直到 13 为止都非常准确，除了 1 之外，所有奇数都应该是 0，但超过 14 的值会给出奇怪的值。例如，当 14 应该给出大约 7/6 的值时，14 给出了一个像 0.9 这样的数字，并且说 22 给出了一个非常负数，大约为 10^-4。奇数给出奇怪的值，比如 15 给出大约 -11。

这里是所有相关代码

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));
}

提前谢谢你。

Answer 1

这里的问题是浮点运算不需要溢出来经历灾难性的精度损失。

浮点数有一个尾数和一个指数，其中数字的值是尾数 * 10^exponent（真正的浮点数使用二进制，我使用十进制）。尾数的精度有限。

当我们添加不同符号的浮点数时，我们最终会得到一个失去精度的最终结果。

例如假设尾数是4位数。 如果我们添加：

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

我们期望得到 1.001 x 10^3。 但是 1.001 x 10^3 + 1.000 x 10^4 = 11.001 x 10^3，表示为 1.100 x 10^4，因为我们的尾数只有 4 位。

所以当我们减去 1.000 x 10^4 时，我们得到 0.100 x 10^4，表示为 1.000 x 10^3 而不是 1.001 x 10^3。

这是一个使用 BigDecimal 的实现，它提供了更好的结果（并且速度慢得多）。

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));
        }
    }
}

尝试更改sum2 中传递给除法的比例并观察对输出的影响。