While computers are very good at doing arithmetic quickly, that speed comes with limitations. One such limitation is handling large numbers, especially when many bits of precision are required. While looking into the collatz problem, what first seemed to be a bug, turned out to be something a bit more interesting.

Initial observation:

While computing values related to collatz for x₀=27 using a C++ program, the value 3⁴¹ was needed. A quick check of log₂ 3⁴¹ shows that 3⁴¹ should require 64.983 bits, so it is a little too big to fit into a 64-bit long long.

To get around the limitations of long long, the long double type was considered for its 112-bit significand. Since 3⁴¹ should fit into 65 bits, the 112-bit significand of long double should be plenty large for it. So, when std::powl(3.0, 41.0) returned an even number, it was most perplexing, as all even numbers have two as a factor.

Question:

Why would std::powl(3.0, 41.0) return an even number when 3⁴¹ must be odd?

Procedure:

To verify this unexpected result, the following simple C++ program was written.

#include <iomanip>
#include <iostream>

int main(int argc, char **argv) {
    long double ThreeTo41 = std::powl(3.0, 41.0);
    std::cout << std::setprecision(20)
              << "3^41: " << ThreeTo41 << std::endl;
    return 0;
}

When run, the code will output:

3^41: 36472996377170786404

however, using an arbitrary precision calculator, it can be determined that the correct output should be:

3^41: 36472996377170786403

which differs by one in the final digit.

The unexpected result is due to Extended Precision. While long double should be 128-bits wide (i.e., sizeof(long double)==16), it is commonly implemented using an 80-bit register. So when 3⁴¹ is computed, which should be about 65 binary or 20 decimal digits long, it should easily fit into a long double’s 112-bit significand. Yet the wrong result is reported.  This is due to the 80-bit implementation only having a 64 bit significand.

To better illustrate the cause of this issue, the table below was generated showing the value computed by std::powl, the expected number of bits needed (i.e., log₂ 3^x), and the hexadecimal representation of the result.

As can be seen, as x increases, the number of bits needed increases, until 3⁴⁰ uses the very last bit of the 64 bit significand.  So, by 3⁴¹, the number no longer fits into 64-bits.

To avoid this problem, libgmp had to be used instead of long double.

Conclusion:

While the long double type is 128 bits, and should have a 112 bit significand, it is often implemented using an 80 bit register with only 64 bits of significand.  This can cause unexpected precision issues when computing larger values.