[R] How to handle large numbers?

[Wacek Kusnierczyk]

Stavros Macrakis wrote:
> On Wed, Feb 11, 2009 at 6:20 AM, Robin Hankin <rksh1 at cam.ac.uk> wrote:
>   
>>> library(Brobdingnag)
>>> exp(1000)/(exp(1007)+5)
>>>       
>> [1] NaN
>>
>>     
>>> as.numeric(exp(as.brob(1000))/(exp(as.brob(1007))+5))
>>>       
>> [1] 0.000911882
>>     
>
> Though brob is certainly useful in many cases, it can't substitute for
> thinking about numeric issues (roundoff, cancellation, overflow,
> underflow) in general.
>
> For example:
>
>   
>> x<-40; log(exp(x)+1)-x
>>     
> [1] 0
>   
>> x<-as.brob(40); log(exp(x)+1)-x
>>     
> [1] -exp(-Inf)
>
> The correct answer is about 4e-18.  Perhaps Ryacas or some other tool
> gets this right, 

bc gets it arbitrarily right:

bc -l <<END
scale=1000
l(e(40)+1)-40
END

as does, e.g., mathematica (which appears to be using the *foss* gmp
library):

N[Log[Exp[40]+1]-40, 1000]


> but in general I don't think it's wise to abdicate
> responsibility to one's tools.

interestingly, perl's bignum gets it slowly and actually wrong (or do i
miss something?):

perl -Mbignum=p,-50 -le 'print (log(exp(40)+1)-40'

vQ