CODEWITHSUNDEEP
r
Mahdi Hadi
Mahdi Hadi

Reputation: 402

Large numbers multiplication in R

I want to know how can I calculate large values multiplication in R. R returns Inf!

For example:

6.350218e+277*2.218789e+215 
[1] Inf

Let me clarify the problem more: consider the following code and the results of outFunc function:

library(hypergeo)
poch <-function(a,b) gamma(a+b)/gamma(a)
n<-c(37 , 41 , 4 , 9 , 12 , 13 , 2 , 5 , 23 , 73 , 129 , 22 , 121 )
v<-c(90.2, 199.3, 61, 38, 176.3, 293.6, 318.6, 328.7, 328.1, 313.3, 142.4, 92.9, 95.5)
DF<-data.frame(n,v)


outFunc<-function(k,w,r,lam,a,b) {
  ((((w*lam)^k) * poch(r,k) * poch(a,b) ) * hypergeo(r+k,a+k,a+b+k,-(w*lam)) )/(poch(a+k,b)*factorial(k))

}

and the function returns:

outFunc(DF$n,DF$v,0.2, 1, 3, 1)
 [1] 0.002911330+  0i 0.003047594+  0i 0.029886646+  0i 0.013560599+  0i 0.010160073+  0i
 [6] 0.008928524+  0i 0.040165795+  0i 0.019402318+  0i 0.005336008+  0i 0.001689114+  0i
[11]         Inf+NaNi 0.005577985+  0i         Inf+NaNi

As can be seen above, outFunc returns Inf+NaNi for n values of 129 and 121. I checked the code sections part by part and I find that the returned results of (wlam)^k poch(r,k) for these n values are Inf. I also check my code with equivalent code in Mathematica which everything is OK:

in:  out[indata[[All, 1]], indata[[All, 2]], 0.2, 1, 3, 1]

out: {0.00291133, 0.00304759, 0.0298866, 0.0135606, 0.0101601, 0.00892852, \
      0.0401658, 0.0194023, 0.00533601, 0.00168911, 0.000506457, \
      0.00557798, 0.000365445}

Now please let me know how we can solve this issue as simple as it is in Mathematica. regards.

Upvotes: 1

Views: 750

Answers (3)

nicola
nicola

Reputation: 24480

For first, I'd recommend two useful reads: logarithms and how floating values are handled by a computer. These are pertinent because with some "tricks" you can handle much bigger values than you think. For instance, your definition of the poch function is terrible. This because the fraction can be simplified a lot but a computer will evaluate the numerator first and if it overflows the result will be useless. That's why R provides beside gamma the lgamma function: it just calculates the logarithm of gamma and can handle much bigger values. So, we calculate the log of each factor in your function and then we use exp to restore the intended values. Try this:

#redefine poch properly
poch<-function(a,b) lgamma(a+b) - lgamma(a)
#redefine outFunc
outFunc<-function(k,w,r,lam,a,b) {
   exp((k*(log(w)+log(lam))+ poch(r,k) + poch(a,b) ) +  
   log(hypergeo(r+k,a+k,a+b+k,-(w*lam)))- poch(a+k,b)-lgamma(k+1))
}
#Now we go
outFunc(DF$n,DF$v,0.2, 1, 3, 1)
#[1] 0.0029113299+0i 0.0030475939+0i 0.0298866458+0i 0.0135605995+0i
#[5] 0.0101600732+0i 0.0089285243+0i 0.0401657947+0i 0.0194023182+0i
#[9] 0.0053360084+0i 0.0016891144+0i 0.0005064566+0i 0.0055779850+0i
#[13] 0.0003654449+0i

Upvotes: 5

Yorgos
Yorgos

Reputation: 30445

> library(gmp)
> x<- pow.bigz(6.350218,277)
> y<- pow.bigz(2.218789,215)
> x*y

Big Integer ('bigz') :
[1] 18592826814872791919942226542714580401488894909642693257011204682802122918146288728149155739011270579948954646130492024596687919148494136290260248656581476275790189359808616520170359345612068099238508437236172770752199936303947098513476300142414338199993261924467166943683593371648

Upvotes: 4

Tim Biegeleisen
Tim Biegeleisen

Reputation: 521194

One option you have available in base R, which does not require a special library, is to convert the two numbers to a common base, and then add the exponents together to get the final result:

> x <- log(6.350218e+277, 10)
> x
[1] 277.8028
> y <- log(2.218789e+215, 10)
> y
[1] 215.3461
> x + y
[1] 493.1489

Since 10^x * 10^y = 10^(x+y), your final answer is 10^493.1489

Note that this solution does not allow to actually store numbers which R would normally treat as INF. Hence, in this example, you still cannot compute 10^493, but you can tease out what the product would be.

Upvotes: 6

Related Questions