Primorial and faster binomial.

[bodrato at mail.dm.unipi.it]

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