1

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?

Community
  • 1
  • 1
Paul Rougieux
  • 10,289
  • 4
  • 68
  • 110
  • You are asking "why" in your question, not "how to solve this". – Roland Nov 04 '16 at 09:57
  • The question title was "is there a way", I edited the question again. It turns out that the `gmp` package doesn't display it correctly either: `library(gmp)` followed by `as.bigz(147573952589676412927)` still returns `Big Integer ('bigz') : [1] 147573952589676412928`. – Paul Rougieux Nov 04 '16 at 10:13
  • 1
    Obviously you need to pass a character and not the floating point number: `as.bigz("147573952589676412927")`. If you pass a floating point number you are passing the imprecise representation. – Roland Nov 04 '16 at 10:21
  • Thanks `as.bigz("147573952589676412927")`, `as.bigz("193707721") * as.bigz("761838257287")` and `as.bigz("2")^67-1` now all display the correct number. – Paul Rougieux Nov 04 '16 at 10:26

1 Answers1

2

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 253 + 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 267

(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".)

Bathsheba
  • 231,907
  • 34
  • 361
  • 483
  • Funny that more than a century after it was made by hand, there still is no straightforward way of making this calculation on a computer. But wait, it does work with python. `pow(2L,67)-1` returns `147573952589676412927L`. I guess an R program could be conceived to display it correctly. – Paul Rougieux Nov 04 '16 at 09:37
  • @PaulRougieux: added some waffle to the answer. – Bathsheba Nov 04 '16 at 09:40
  • Thanks for mentioning the Brobdingnag package. [This answer](http://stackoverflow.com/a/32368658/2641825) mentions the [gmp package](https://cran.r-project.org/package=gmp) which seems to have been updated more recently. – Paul Rougieux Nov 04 '16 at 09:46