Karim Belabas on Sun, 27 Nov 2022 10:59:47 +0100


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

Re: binomial() challenge (+thanks) 

* Ruud H.G. van Tol [2022-11-26 13:17]:
> 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))
[...]
> ? 1.0* binomial(10^40,43)
> %16 = 1.6552108677421951886982956337138493958 E1667
> 
> ? 1.0* binom_2(10^40,43)
> %17 = 7.416480782428988905 E1778

This is expected of course: since n is huge and k small, we have n+1 ~ n-k+1
and the subtraction lngamma(n+1) - lngamma(n-k+1) suffers from massive
cancellation.

(10:44) gp > n=10^40;k=43
%1 = 43
(10:44) gp > lngamma(n+1) - lngamma(n-k+1)
%2 = 4096.000000000000000

But in fact
(10:44) gp > \p100
   realprecision = 115 significant digits (100 digits displayed)
(10:44) gp > lngamma(n+1) - lngamma(n-k+1)
%3 = 3960.446359949758576510945302[...]

So binom_2() is a terrible idea without complicated counter-measures. On the
other hand, binom_1() is fine since divisions do not cause cancellation.
Actually, gamma(n+1) / (gamma(k+1) * gamma(n-k+1)) would be just as stable
and marginally faster (you may use factorial(), as Bill suggested).

> So it looks to me, that always using binomial(), so never even
> consider using an alternative, is simply the right way.

Not quite. :-)

binomial(n,k) is exact, so very wasteful if only a floating point
approximation is required. The underlying algorithm (when n fits into an
unsigned long) actually computes the factorisation of the result using
Lucas's formula for the valuation at each prime. This is then expanded
as an exact (and possibly huge!) integer.

It's not straightforward for a user to re-implement this in GP but it's
easy to adapt the C code so as to return the binomial in factored form
or expand it as a floating point number. Both would be massively more
efficient than binomial(). Lacking applications, I won't do it, but I'm
dropping the comment here in case anybody is interested.

Cheers,

    K.B.
--
Karim Belabas  /  U. Bordeaux,  vice-président en charge du Numérique
Institut de Mathématiques de Bordeaux UMR 5251 - (+33) 05 40 00 29 77
http://www.math.u-bordeaux.fr/~kbelabas/
`