Binomial improvements
Torbjorn Granlund
tg at gmplib.org
Sun Dec 11 21:15:29 CET 2011
nisse at lysator.liu.se (Niels Möller) writes:
At least a table för k should take little memory. For each k < k', it
could tabulate not only k! but also a hensel inverse (of the odd part,
and the shift count) of one or two limbs. Should be useful also for the
case of large n and small k.
I implemented this, see https://gmplib.org/devel/.
The function becomes extremely simple:
void
mpz_bc_bin_uiui (mpz_ptr r, unsigned long int n, unsigned long int k)
{
mp_limb_t i, j, c;
c = n - k + 1;
i = n - k + 2;
for (j = 1; j < k; j++)
{
c = c * i;
i++;
c = c * zinv[j] >> ctz[j];
}
PTR(r)[0] = c;
SIZ(r) = 1;
}
There are also tables, currently 32 limbs large, tabulating values also
for even numbers. That might be a bit silly. One might as well
tabulate just odd numbers, and ctz[j] to the index, like this
(untested):
for (j = 1; j < k; j++)
{
c = c * i;
i++;
c = c * zinv[j >> ctz[j]] >> ctz[j];
}
Is this the best way to handle really small numbers? This code cannot
handle k > 32 (with present tables) and n greater then GMP_LIMB_BITS +
2.
One may write a variant that uses mpn_mul_1 instead, to allow any n that
fits in a limb.
Unlike my bdiv_bin_uiui code, this code puts all multiplies and
"divides" on the same dependency chain. This might hurt performance on
some machines (with good multiply throughput and long latency).
Clearly, pre-generating the actual results for small n,k is faster....
--
Torbjörn
More information about the gmp-devel
mailing list