In R why is factorial(100) displayed differently to prod(1:100)?

Question

In R I am finding some odd behaviour that I can't explain and I am hoping someone here can. I believe that the value of 100! is this big number.

A few lines from the console showing expected behaviour...

>factorial( 10 ) [1] 3628800 >prod( 1:10 ) [1] 3628800 > prod( as.double(1:10) ) [1] 3628800 > cumprod( 1:10 ) [1]       1       2       6      24     120     720    5040   40320  362880 3628800 

However when I try 100! I get (notice how the resulting numbers begin to differ about 14 digits in):

> options(scipen=200) #set so the whole number shows in the output > factorial(100) [1] 93326215443942248650123855988187884417589065162466533279019703073787172439798159584162769794613566466294295348586598751018383869128892469242002299597101203456 > prod(1:100) [1] 93326215443944102188325606108575267240944254854960571509166910400407995064242937148632694030450512898042989296944474898258737204311236641477561877016501813248 > prod( as.double(1:100) ) [1] 93326215443944150965646704795953882578400970373184098831012889540582227238570431295066113089288327277825849664006524270554535976289719382852181865895959724032 > all.equal( prod(1:100) , factorial(100) , prod( as.double(1:100) ) ) [1] TRUE 

If I do some tests against a variable set to the 'known' number of 100! then I see the following:

# This is (as far as I know) the 'true' value of 100! > n<- as.double(93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000) > factorial(100) - n [1] -1902315522848807765998160811905210717565551993186466795054798772271710903343294674760811531554315419925519536152107160826913610179566298858520576 > prod(1:100) - n [1] -48777321098687378615337456715518223527321845979140174232174327494146433419058837814379782860367062049372295798771978482741374619988879457910784 > prod(as.double(1:100)) - n [1] 0 

The final result evaluates to zero, but the number returned for prod( as.double( 1:100 ) ) does not display as I would expect, even though it correctly evaluates prod( as.double( 1:100 ) ) - n where n is a variable set to the value of 100!.

Can anyone explain this behaviour to me please? It should not be related to overflow etc as far as I am aware, as I am using a x64 system. Version and machine info below:

> .Machine$double.xmax [1] 1.798e+308 > str( R.Version() ) List of 14  $ platform      : chr "x86_64-apple-darwin9.8.0"  $ arch          : chr "x86_64"  $ os            : chr "darwin9.8.0"  $ system        : chr "x86_64, darwin9.8.0"  $ status        : chr ""  $ major         : chr "2"  $ minor         : chr "15.2"  $ year          : chr "2012"  $ month         : chr "10"  $ day           : chr "26"  $ svn rev       : chr "61015"  $ language      : chr "R"  $ version.string: chr "R version 2.15.2 (2012-10-26)"  $ nickname      : chr "Trick or Treat" 

Can anyone explain this to me? I don't doubt that R does everything correctly and this is most likely useR related. You might point out that since prod( as.double( 1:100 ) ) - n evaluates correctly what am I bothered about, but I am doing Project Euler Problem 20 so I needed the correct digits displayed.

Thanks

Answer 1

This has to do not with the maximum value for a double but with its precision.

100! has 158 significant (decimal) digits. IEEE doubles (64 bit) have 52 bits of storage space for the mantissa, so you get rounding errors after about 16 decimal digits of precision have been exceeded.

Incidentally, 100! is in fact, as you suspected,

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000 

so all of the values R calculated are incorrect.

Now I don't know R, but it seems that all.equal() converts all three of those values to floats before comparing, and so their differences are lost.