CODEWITHSUNDEEP
rfactorial
Valegard234
Valegard234

Reputation: 45

Factorial(x) for x>170 using Rmpfr/gmp library

The problem that I would like to solve is the infinite sum over the following function:

enter image description here

For the sum I use an FTOL determination criterion. This whole term doesn't create any problems until z becomes very large. I expect the maximum value of z around 220. As you can see the first term has its max around Factorial(221) and therefore has to go around Factorial(500) until the determination criterion has been reached. After spotting this problem I didn't want to change the whole code (as it is only one small part) and tried to use library('Rmpfr') and library('gmp'). The problem is that I do not get what I want to. While multiplication normally works, subtraction fails for higher values:

This works

> factorialZ(22)-factorial(22)
Big Integer ('bigz') :
[1] 0

but this fails: 

> factorialZ(50)-factorial(50)
Big Integer ('bigz') :
[1] 359073645150499628823711419759505502520867983196160

another way I tried: 

> gamma(as(10,"mpfr"))-factorial(9)
1 'mpfr' number of precision  128   bits 
[1] 0

> gamma(as(40,"mpfr"))-factorial(39)
1 'mpfr' number of precision  128   bits 
[1] 1770811808798664813196481658880

There has to be something that I don't really understand. Does someone have a even better solution for the problem or can someone help me out with the issue above?

Upvotes: 1

Views: 460

Answers (2)

oropendola
oropendola

Reputation: 1101

50! is between 2^214 and 2^215 so the closest representable numbers are 2^(214-52) apart. factorial in R is based on a Lanczos approximation whereas factorialZ is calculating it exactly. The answers are within the machine precision:

> all.equal(as.numeric(factorialZ(50)), factorial(50))
[1] TRUE

The part that you're not understanding is floating point and it's limitations. You're only getting ~15 digits of precision in floating point. factorialZ(50) has a LOT more precision than that, so you shouldn't expect them to be the same.

Upvotes: 2

Carl Witthoft
Carl Witthoft

Reputation: 21502

I think you incorrectly understand the priorities in factorialZ(x)-factorial(x) . The second term, factorial(x) is calculated before it's converted to a bigz to be combined with the first term.

You must create any integer outside the 2^64 (or whatever, depending on your machine) range using a bigz - compatible function.

Upvotes: 2

Related Questions