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

Question

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

在 13 之前都非常准确，除 1 外所有奇数都应为 0，但超过 14 的值会给出奇怪的值。例如，14 给出一个像 0.9 这样的数字，而它应该给出大约 7/6 的值，而 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^指数（真正的浮点数使用二进制，我使用十进制）。尾数精度有限。

当我们添加不同符号的浮点数时，我们最终的结果可能会丢失精度。

例如假设尾数是 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 并观察对输出的影响。

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

Issues with accuracy with Bernoulli Number generator in java

java

math

for-loop

sum

bernoulli-numbers