Почему моя программа для вычисления чисел Вудолла дает неверные результаты после n> 47?

Question

Для этой функции, которая вычисляет числа Вудолла до n = 64

А алгоритм для дерева: W _n = n ⋅ 2 ⁿ - 1

for (int n = 1; n <= 64; ++n)
{
    a[n - 1] = (n * (exp2(n))) - 1;
}

Но после того, как n превысит 47, результаты окажутся неверными в том смысле, что кажется, что он забывает - 1 результат n * (exp2(n)).

Вот что выводится, если я cout значения через

std::cout << i << ":\t" << std::setprecision(32) << a[i - 1] << std::endl;

... прежде чем правильно

n
45:     1583296743997439
46:     3236962232172543
47:     6614661952700415
48:     13510798882111488
49:     27584547717644288
50:     56294995342131200

... после неверного

для a[] является длинным без знака int

Функция выдаст правильные результаты, если я отделю операцию - 1 от ее собственного цикла for:

for (int n = 1; n <= 64; ++n)
{
    a[n - 1] = (n * (exp2(n)));
}


for (int n = 1; n <= 64; ++n)
{
    a[n - 1] = a[n - 1] - 1;
}

Answer 1

double имеет ограничения по точности. Он использует двоичную базу для работы, хотя, означает, что большинство чисел, заканчивающихся серией нулевых битов в двоичном коде, может быть представлено точно, что имеет место для кратных exp2(int).

50 * exp2(50), например, 56294995342131200, это C8000000000000 в шестнадцатеричном формате. Даже если количество цифр превышает ограничения точности double, оно может быть представлено точно. Однако, если я попытаюсь сложить или вычесть 1 из этого числа, это уже не так.

double не может представлять 56294995342131199 или 56294995342131201, поэтому при попытке сделать это просто округляется до 56294995342131200.

Вот почему ваш бит - 1 не работает, он все еще работает как double, когда вы пытаетесь выполнить эту операцию. Вам придется привести остаток выражения к int64_t перед выполнением этого вычитания.

Но другое решение - вообще не использовать exp2(). Поскольку мы работаем с целыми числами, вы можете просто использовать побитовые операции для выполнения той же задачи. (1 << n) даст вам те же результаты, что и exp2(), за исключением того, что он теперь в целочисленном формате, и, поскольку вы просто умножаете это значение на n, вы можете просто сделать (n << n).

Конечно, это все равно сломает черту. int64_t может содержать только такие большие числа, как 2 ⁶³ -1 и uint64_t 2 ⁶⁴ -1, которые должны прерваться, когда ваш итератор достигнет n = 57.

Answer 2

exp2(n) возвращает double.

В IEEE754 (очень распространенная спецификация для типов с плавающей запятой), которая дает вам только точные целые числа до 52-й степени 2. После этого вы получите приближения.

Вы наблюдаете проблемы до 52-го числа Вудолла, поскольку все выражение n * (exp2(n))) - 1 равно double из-за неявного преобразования типов. По вычислительной причуде, именно -1 вызывает проблему. Просто случается, что другой термин является подходящим кратным степени 2, что позволяет представить его как двойное без потери точности! Это причина того, что ваш второй фрагмент работает, а ваш первый - нет.

В системе с 64-битным int вы достигнете целочисленных пределов (и неопределенного поведения) при 63-й степени 2 *.

Лучше всего генерировать числа Вудолла чисто в unsigned арифметике (обратите внимание на соотношение между << и степенью 2), возможно, даже используя рекуррентное отношение для последовательных чисел Вудолла.