Binomial -- faster for large args

[abbott at module.dima.unige.it]

Hi,

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 -- I didn't compare then for
smaller values.

I have NOT put in special case code (e.g. r=0 or r=1 or r > n/2)

Hope it is clean enough for you to use.
John.
---------------------------------------------------------------------
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_mul(ans,ans,num);
    mpz_divexact(ans,ans,den);
  }
}