Binomial -- faster for large args
abbott at module.dima.unige.it
abbott at module.dima.unige.it
Mon May 17 13:13:06 CEST 2004
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);
}
More information about the gmp-devel
mailing list