[R] About 'choose' function

[(Ted Harding)]

On 08-Nov-04 Ted Harding wrote:
> On 08-Nov-04 Ted Harding wrote:
>> [...]
>> 
>>   > choose(60000000000,31)
>>   [1] 1.612899e+300
>> 
>> [...]
>> 
>>   my.choose3<-function(x,y){
>>     if((x==y)||(y==0)) return(1);
>>     m <- min(y,x-y)
>>    prod(seq(x,(x-m+1),by=-1)/(seq(m,1,by=-1)))
>>   }
>> 
>> when
>> 
>>   my.choose3(60000000000,31)
>>   [1] 1.613121e+300
>> 
>> which has a slight difference (0.014% greater) from the result of
>> 'choose'.
>> 
>> Given the method of computation, I might feel inclined to trust
>> 'my.choose3' rather than 'choose', but I'm not sure of this
>> without studying the internal code of 'choose', and would welcome
>> comments!
> 
> I just did this calculation using 'bc' (which works to arbitrary
> decimal precision), getting (after rounding 59999999999.99...996,
> 1744 9's in all, to 60000000000):
> 
> 16131211663389665874291103209859887097643782097733433535176164529013\
> 93347758720756812586746724689691665340326913939172715147819303293449\
> 54257506586926613288215108141984988883561331127046497391143349652601\
> 81932350377398578090978006727056282867666693363201127636940304448887\
> 82022049197129032260000000000
> 
> so the leading digits are 1613121..., agreeing with 'my.choose3'
> as above, rather than with 'choose'!

Further investigation indicates that there is an integer truncation
problem with 'choose':

  > print(choose(60000000000,1),digits=20)
  [1] 60001679906.00000
  > print(my.choose3(60000000000,1),digits=20)
  [1] 6e+10

Best wishes to all,
Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861  [NB: New number!]
Date: 08-Nov-04                                       Time: 15:56:49
------------------------------ XFMail ------------------------------