Primorial and faster binomial.
bodrato at mail.dm.unipi.it
bodrato at mail.dm.unipi.it
Thu Apr 19 12:10:46 CEST 2012
Ciao!
I pushed some new code in the main repository:
- a prime-factor-based implementation for the binomial,
- a new function to compute the primorial.
Both are based on the bit-wise prime-sieve I wrote for the factorial,
that's why I moved it in the main directory, it is now called
gmp_primesieve, but this name might be a bad idea, since it does not use
the gmp_primesieve_t type...
The current mpz_primorial_ui function is almost trivial, it sieves all
primes up to the given limit (storing the full bit-array), then it
produces a list of limb-sized factors and calls mpz_prodlimbs. It
implements one of the possible definitions of primorial(N): product of
primes <= N, NOT the product of the N smaller primes...
The new code for binomial(N,K) should be fast asymptotically, i.e. it
should give very good results with tune/speed mpz_bin_uiui . But it
obviously is not always fast: as usual we need a threshold mechanism, but
we are in a bi-dimensional problem.
A threshold on K is appealing, but incorrect. The prime-factor-based
implementations sieves primes up to N, i.e. with a fixed K its complexity
grows at least linearly with N. The complexity of the basecase algorithm
grows with log^2(N)...
The current code uses the basecase below an hard-coded threshold on K, and
anytime K < N/16. Tuning, and a deeper analysis are needed.
There still are a lot of duplicated tables, macros and small functions,
that should be shared among fac_ui, bin_uiui, primorial_ui... Two, in
particular, may be interesting in general (or duplicated!): limb_apprsqrt,
and
/* Returns an approximation of the sqare root of x. *
* It gives: x <= limb_apprsqrt (x) ^ 2 < x * 9/4 */
static mp_limb_t
limb_apprsqrt (mp_limb_t x)
{
int s;
ASSERT (x > 2);
count_leading_zeros (s, x - 1);
s = GMP_LIMB_BITS - 1 - s;
return (CNST_LIMB(1) << (s >> 1)) + (CNST_LIMB(1) << ((s - 1) >> 1));
}
/* n^log <= GMP_NUMB_MAX, a limb can store log factors less than n */
static unsigned
log_n_max (mp_limb_t n)
{
static const mp_limb_t table[] = { NTH_ROOT_NUMB_MASK_TABLE };
unsigned log;
for (log = numberof (table); n > table[log - 1]; log--);
return log;
} /* NTH_ROOT_NUMB_MASK_TABLE is generated by gen-fac_ui and stored in
fac_ui.h */
I'd appreciate comments on any of the recent changes.
Best regards,
Marco
--
http://bodrato.it/software/combinatorics.html
More information about the gmp-devel
mailing list