Search code examples
rmathfactorial

R factorial() is not yielding the correct result


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?


Solution

  • 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