Ruud H.G. van Tol on Sat, 26 Nov 2022 14:56:05 +0100


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: binomial() challenge (+thanks) 


On 2022-11-26 14:29, Gottfried Helms wrote:

Am 26.11.2022 um 13:17 schrieb Ruud H.G. van Tol:



binom_2(n,k) = exp(lngamma(n+1) - lngamma(k+1) - lngamma(n-k+1))

? version()
%12 = [2, 15, 1]

? 1.0* binomial(10^40,43)
%16 = 1.6552108677421951886982956337138493958 E1667

? 1.0* binom_2(10^40,43)
%17 = 7.416480782428988905 E1778


Hmm, this surprised me, because I've never observed
such discrepancies, although I've experimented with
such implementations.
However, I work by default with dec precision of
200 digits, and reproducing your function calls
led to correct digits as well for the binom_2() version:

  \\ -----------------------------------------------------------------------------------------------
  version() \\ %75 = [2, 16, 0, 28083, "5ad72a2e8c"]

  fmt(200,60) \\ "prec=211 display=0"

  1.0* binomial(2^55,43) \\ %101 = 1.42821330463476670428642354475995521176331125070436530629880 E659
  1.0* binom_1(2^55,43)  \\ %103 = 1.42821330463476670428642354475995521176331125070436530629880 E659
  1.0* binom_2(2^55,43)  \\ %105 = 1.42821330463476670428642354475995521176331125070436530629880 E659

  1.0* binomial(10^40,43) \\ %107 = 1.65521086774219518869829563371384939580110837820413175684538 E1667
  1.0* binom_2(10^40,43)  \\ %109 = 1.65521086774219518869829563371384939580110837820413175684538 E1667

  \\ ------------------------------------------------------------------------------------------------
(even to 120 dec digits no difference)
Thanks. An implied message was, that binomial itself is in general the best one to use. Also try 2^62 and 10^60 and such. 

-- Ruud