2

I am trying to calculate the factorial of 52 in R. Astonishingly, I am getting contradicting results.

aaa<-1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20*21*22*23*24*25*26*27*28*29*30*31*32*33*34*35*36*37*38*39*40*41*42*43*44*45*46*47*48*49*50*51*52

bbb<-factorial(52)

aaa
[1] 80658175170943876845634591553351679477960544579306048386139594686464

bbb
[1] 80658175170944942408940349866698506766127860028660283290685487972352

aaa==bbb #False

What am I doing wrong?

enroute
  • 183
  • 9
  • To make my confusion complete: This site gives another solution (80658175170943878571660636856403766975289505440883277824000000000000) https://czep.net/weblog/52cards.html – enroute Feb 17 '23 at 11:55
  • I guess it's a floating point precision issue. – harre Feb 17 '23 at 12:04
  • 3
    [Why are these numbers not equal?](https://stackoverflow.com/q/9508518/10068985) – Darren Tsai Feb 17 '23 at 12:11
  • 3
    In R, the maximum integer size is 2^31 - 1. Therefore any factorial larger than `factorial(12)` is too large to be stored as an integer and has to be stored as a floating point number. Even with 64-bit systems, this means you will only have at most 18 significant digits. When you multiply floating point numbers together, this inaccuracy is compounded, so different ways of calculating the same number will give slightly different results. The two numbers here are equal within 1 part in 10 trillion, which really isn't too bad. – Allan Cameron Feb 17 '23 at 12:12

1 Answers1

4

This is a well known problem in computing with large numbers; R uses double-precision floating-point, the precision of which may vary by machine. Thats why you are getting multiple results across methods (including the online calculator in your comments). If you want to change your precision (in bits), one option is to use the Rmpfr package:

Rmpfr::mpfr(factorial(52), 
          6) # six bits

#1 'mpfr' number of precision  6   bits 
#[1] 8.09e+67

Rmpfr::mpfr(factorial(52), 
          8) # eight bits

#1 'mpfr' number of precision  8   bits 
#[1] 8.046e+67

This will allow you to obtain a result with the same value:

x <-Rmpfr::mpfr(1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20*21*22*23*24*25*26*27*28*29*30*31*32*33*34*35*36*37*38*39*40*41*42*43*44*45*46*47*48*49*50*51*52, 
          8)
y <- Rmpfr::mpfr(factorial(52), 
                 8)

x == y
#[1] TRUE
jpsmith
  • 11,023
  • 5
  • 15
  • 36