Why is my program for calculating the woodall numbers producing wrong results after n >47?

Question

For this function that calculates the woodall numbers up to n = 64

And the algorithm for a woodall is W_n = n ⋅ 2ⁿ - 1

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

But after n is greater than 47 the results are wrong in that it seems like it is forgetting to - 1 the result of n * (exp2(n)).

Here is what the output is if i cout the values via

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

... before is correct

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

... after is incorrect

for a[] is an unsigned long int

The function produces correct results if I separate the - 1 operation out to its own for loop though:

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

exp2(n) returns a double.

In IEEE754 (a very common specification for floating point types), that only gives you exact integers up to the 52nd power of 2. Thereafter you get approximations.

You observe issues before the 52nd Woodall number since the entire expression n * (exp2(n))) - 1 is a double due to implicit type conversion. By a computational quirk, it's the -1 that causes the problem. It just happens that the other term is an appropriate multiple of a power of 2 which allows it to be represented as a double without precision loss! This is the reason behind your second snippet working but your first snippet not.

On a system with a 64 bit int, you'll hit integer limits (and undefined behaviour) on the 63rd power of 2.

Your best bet is to generate the Woodall numbers purely in unsigned arithmetic (note the relationship between << and a power of 2), perhaps even using a recurrence relation for successive Woodall numbers.

Answer 2

double has precision limitations. It does use a binary base to work, though, meaning most numbers finishing with a series of zero bits in binary can be represented exactly, which is the case for multiples of exp2(int).

50 * exp2(50) which is 56294995342131200 for example, is C8000000000000 in hexadecimal. Even though the number of digits exceeds the precision limitations of double, it can be represented exactly. However, if I try to sum or subtract 1 from this number, that is no longer the case.

double can't represent 56294995342131199 nor 56294995342131201, so when you try to do it, it simply gets rounded back to 56294995342131200.

This is why your - 1 bit is failing, it is still being operated as a double when you try to perform this operation. You'd have to cast the rest of the expression to int64_t before performing this subtraction.

But another solution is to not use exp2() at all. Since we are working with integers, you can simply use bitwise operations to perform the same task. (1 << n) will yield you the same results as exp2() except it is now in integer format, and because you are just multiplying this to n, you can actually just do (n << n).

Of course, this will still break down the line. int64_t can only hold a number as big as 2⁶³-1 and uint64_t 2⁶⁴-1, which should break when you iterator reaches around n = 57.