2.1 Factorial Function

Truncation error and round-off error can also be occurred in implementing the factorial function.

Sterling’s formula for large value approximation is quite inaccurate, and its use should be discouraged now. Lanczos approximation is much better to use now practically.

naiveFactorial.a = function(x) # do not use for real work situation
{ # x should be a positive integer.
  Result = 1
  for (i in 1:x) Result = Result*i
  return(Result)
}
format(naiveFactorial.a(170), digits=22)

## [1] "7.2574156153079941358e+306"

format(prod(1:170), digits=22)

## [1] "7.257415615307999125e+306"

naiveFactorial.a(170) == prod(1:170)

## [1] FALSE

identical(naiveFactorial.a(170), prod(1:170)) # bitwise comparison

## [1] FALSE

all.equal(naiveFactorial.a(170), prod(1:170)) # practically same?

## [1] TRUE

naiveFactorial.b = function(x) # Do not use for real work situation
{ # x should be a positive integer.
  Result = 1
  for (i in x:1) Result = Result*i
  return(Result)
}
format(naiveFactorial.b(170), digits=22)

## [1] "7.257415615308004115235e+306"

format(prod(170:1), digits=22)

## [1] "7.257415615307999125e+306"

prod(1:170) == prod(170:1) # TRUE in R 3.5.1

## [1] TRUE

betterFactorial.a = function(x)
{ # x should be a positive integer.
  Result = prod(seq(1, x, 2))*prod(seq(2, x, 2))
  return(Result)
}
format(betterFactorial.a(170), digits=22)

## [1] "7.257415615307999125e+306"

prod(1:170) == betterFactorial.a(170)

## [1] TRUE

Practically fast and precise one.

format(tableFactorial(170), digits=22)

## [1] "7.257415615307999125e+306"

Comparison of testFactorial function and factorial R internal function.

Start = 1
End = 170      # maximum value not to be infinity
RelTol = 1e-13 # relative tolerance
for (i in Start:End) {
  x1 = tableFactorial(i)
  x2 = factorial(i)
  RelDiff = abs((x1 - x2)/x2) # relative difference
  if (RelDiff > RelTol) print(paste(i, format(c(x1, x2), digits=22)))
}

## [1] "127 3.012660018457659457086e+213" "127 3.012660018457322778484e+213"
## [1] "132 1.118248651196004340015e+224" "132 1.118248651196143034952e+224"
## [1] "134 1.992942746161518811215e+228" "134 1.992942746161724840246e+228"
## [1] "151 8.627209774233240077712e+264" "151 8.627209774234339096158e+264"
## [1] "153 2.006343905095682317264e+269" "153 2.006343905095905783042e+269"
## [1] "154 3.089769613847350809366e+271" "154 3.089769613847712671816e+271"
## [1] "166 9.003691705778437654600e+297" "166 9.003691705779608680550e+297"
## [1] "170 7.257415615307999125303e+306" "170 7.257415615308882284734e+306"

If you make RelTol more strict, you can see more values. Values from tableFactorial are more precise. In fact, results of R internal function prod are same to tableFactorial function.

for (i in 2:170) if (prod(1:i) != tableFactorial(i)) print(i) # nothing to print