Проблемы с точностью с помощью генератора чисел Бернулли в java

Question

Я создал код, который генерирует числа Бернулли на основе формулы 33 в MathWorld. Это приведено на https://mathworld.wolfram.com/BernoulliNumber.html и должно работать для всех целых чисел n, но оно очень быстро расходится с ожидаемыми результатами, как только достигает n= 14. Я думаю, что проблема может быть в коде факториала, хотя я понятия не имею.

Это довольно точно до 13, все нечетные числа должны быть 0, кроме 1, но значения после 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 х 10^3 + 1.000 х 10^4 - 1.000 х 10^4

мы ожидаем получить 1,001 x 10^3. Но 1,001 х 10 ^ 3 + 1,000 х 10 ^ 4 = 11,001 х 10 ^ 3, что представляется как 1,100 х 10 ^ 4, учитывая, что наша мантисса имеет только 4 цифры.

Поэтому, когда мы вычитаем 1,000 x 10^ 4, мы получаем 0,100 x 10^ 4, что представляется как 1,000 x 10^ 3, а не 1,001 x 10^3.

Вот реализация с использованием что дает лучшие результаты (и намного медленнее).

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

Попробуйте изменить масштаб, переданный делению в и наблюдайте за эффектом на выходе.