5

I was showing my daughter the magic pattern of two 111..111 number multiplication in R, e.g. 

> options(scipen=999)
> 1^2
[1] 1
> 11^2
[1] 121
> 111^2
[1] 12321
> 1111^2
[1] 1234321
> 11111^2
[1] 123454321
> 111111^2
[1] 12345654321
> 1111111^2
[1] 1234567654321
> 11111111^2
[1] 123456787654321
> 111111111^2
[1] 12345678987654320

Everything works perfectly until it goes to 9 digits. As you can tell from the last line, the answer is obviously wrong; it should be 1234567898765432*1*, not 1234567898765432*0*.

I am running this in R (v3.6.1, x86_64-apple-darwin15.6.0, 64-bit) on a Macbook Pro laptop.

Someone said that it might be caused by an integer overflow. Two questions here:

  1. Can anyone explain how this happened exactly? e.g. how was the last 1 changed to 0?

  2. How can I get the correct answer in R?

Thanks.

Xianjun
  • 161
  • 6
  • 3
    try using the `bit64` package for large integer representation: https://cran.r-project.org/web/packages/bit64/index.html – mikeck Oct 20 '19 at 19:20
  • 1
    you have to do `bit64::as.integer64('111111111') * bit64::as.integer64('111111111')`, though; `^.integer64` coerces to double – alistaire Oct 20 '19 at 19:28
  • 1
    or gmp: `gmp::as.bigz('111111111')^2`. Relevant: https://stackoverflow.com/questions/2053397/long-bigint-decimal-equivalent-datatype-in-r – alistaire Oct 20 '19 at 19:36
  • Thanks to @alistaire for providing the solution to getting the correct answer. Can anyone explain how exactly the last 1 was changed to 0 in the normal R? – Xianjun Oct 21 '19 at 01:32
  • If you call `str` on the result, even if you start with integers `str(111111111L^2L)`, you end up with a double, so it's floating point error. If you use multiplication instead of a power with integers `111111111L * 111111111L`, you get integer overflow, which is a warning; it returns `NA`. – alistaire Oct 21 '19 at 03:45

1 Answers1

2

My personal favorite for this kind of problem is package Rmpfr, an R package for the GNU MPFR Library.

library(Rmpfr)

## 53 bit precision
ones9 <- mpfr(1111111111, precBits = 53)

ones9^2
#1 'mpfr' number of precision  53   bits 
#[1] 1234567900987654400


## 100 bit precision
ones9b <- mpfr(1111111111, precBits = 100)
ones9b^2
#1 'mpfr' number of precision  100   bits 
#[1] 1234567900987654321
Rui Barradas
  • 70,273
  • 8
  • 34
  • 66