Gottfried Helms on Sat, 26 Nov 2022 14:31:00 +0100


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

Re: binomial() challenge (+thanks)


Am 26.11.2022 um 13:17 schrieb Ruud H.G. van Tol:
> 
> PARI binomial() challenge:
> 
> binom_1(n,k) = gamma(n+1) / gamma(k+1) / gamma(n-k+1)
> binom_2(n,k) = exp(lngamma(n+1) - lngamma(k+1) - lngamma(n-k+1))
> 
> ? version()
> %12 = [2, 15, 1]
> 
> ? 1.0* binomial(2^55,43)
> %13 = 1.4282133046347667042864235447599552118 E659
> 
> ? 1.0* binom_1(2^55,43)
> %14 = 1.4282133046347667042864235447599552118 E659
> 
> ? 1.0* binom_2(2^55,43)
> %15 = 1.4282133046347667042944763056801257481 E659
> 
> and then
> 
> ? 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)

Gottfried Helms