Paul Rougieux - 2 months ago 10

R Question

In 1876, Frank Nelson Cole showed that the Mersenne number M67 = 2^67 − 1 is not a prime number.

During Cole's so-called "lecture", he approached the chalkboard and in complete silence proceeded to calculate the value of M67, with the result being 147,573,952,589,676,412,927. Cole then moved to the other side of the board and wrote 193,707,721 × 761,838,257,287, and worked through the tedious calculations by hand. Upon completing the multiplication and demonstrating that the result equaled M67, Cole returned to his seat, not having uttered a word during the hour-long presentation. His audience greeted the presentation with a standing ovation.

I tried to calculate this with R. But it seems there is a rounding error. (This answer explains how to show more digits.)

`print(2^67-1, digits =21)`

# [1] 147573952589676412928

print(193707721 * 761838257287, digits = 21)

# [1] 147573952589676412928

print(147573952589676412927, digits = 21)

# [1] 147573952589676412928

Is there a way to prevent R from rounding at the next integer?

Answer

You've hit a limit of 64 bit double precision floating point.

IEEE754 double precision floating point can only display integers exactly up to the 53rd power of 2. Thereafter it will round to the nearest available integer. You can see this for yourself by considering the first number that exhibits this:

`9,007,199,254,740,993`

is 2^{53} + 1 and will be returned as `9,007,199,254,740,992`

Your number, `147,573,952,589,676,412,927`

, behaves similarly and will snap to the adjacent 2^{67}

(You need to use a library that can handle big integers. Java has a good one although the lack of operator overloading in Java makes it a pig to use. Many C++ compilers have a 128 bit integer type which will be good enough for this. I think even R has a large number library - Brobdingnag? - named, somewhat pretentiously, after one of the islands in Jonathan Swift's "Gulliver's Travels".)