Binomial -- faster for large args

[abbott at module.dima.unige.it]

On Sat, 15 May 2004, Kevin Ryde wrote:
> abbott at module.dima.unige.it writes:
> > The following routine (called "binomial") seems to be usefully
> > faster than that supplied with GMP 4.1.3 (on my Linux/Athlon box)
> > at least for larger values of n and r
>
> gcd is normally pretty slow, without actually measuring it I'd imagine
> just dividing would be better.
>
> Oh, and divexact is not yet sub-quadratic, plain mpz_tdiv_q might be
> faster.

I would expect the divisions to be quite unbalanced most of the time,
so it may not matter much that divexact is "naive".  I did try a more
balanced version but it didn't work as well -- maybe there was a bug
(or it could have been the "naive" nature of divexact?).

A quick check here shows that the cross-over is pretty high:
my routine is faster if the larger arg to "binomial" is 5000 or more.
For still larger arguments the gain is more noticeable:
  e.g. to compute binomial(1000000, 500000)
       mpz_bin_uiui  takes about       100s CPU
       binomial (my code) takes about   13s CPU

Maybe no one will ever need to compute such large values; still, the code
might be useful as a benchmark for future improvements to mpz_bin_uiui.

John.
PS below is corrected version of the code I sent last week
-------------------------------------------------------------
void binomial(mpz_t ans, unsigned long n, unsigned long r)
{
  /* The following value for chunk worked well in my tests */
  const unsigned chunk = 10+int(floor(2*sqrt(double(n))));

  mpz_t num, den, g;
  mpz_init(num);
  mpz_init(den);
  mpz_init(g);
  mpz_set_ui(ans,1);
  unsigned int k=0;
  while (k < r)
  {
    mpz_set_ui(num,1);
    mpz_set_ui(den,1);
    for (unsigned int i=0; i < chunk; ++i)
    {
      mpz_mul_ui(num,num,n-k);
      mpz_mul_ui(den,den,k+1);
      ++k;
      if (k >= r) break;
    }
    mpz_gcd(g,num,den);
    mpz_divexact(num,num,g);
    mpz_divexact(den,den,g);
    mpz_divexact(ans,ans,den);
    mpz_mul(ans,ans,num);
  }
  mpz_clear(g);
  mpz_clear(den);
  mpz_clear(num);
}