Double factorial and primorial

[Torbjorn Granlund]

nisse at lysator.liu.se (Niels Möller) writes:

  > In fac_ui.c there already is a function for sieving: bitwise_primesieve.
  > It fills an mp_ptr "with the characteristic function of composite"s, and
  > it returns the count of "prime integers in the interval [4, n]".
  
  That's going to use a lot more memory than necessary. To sieve all
  primes up to n, you need to first create a list of primes up to sqrt(n),
  and represent it in any convenient way which lets you iterate over it
  (bitmap, or vector of differences, being the two most obvious choices).
  
  Above sqrt(n), do the sieving blockwise. Create a bitmap of a suitable
  blocksize (it should fit in the cache), let this bitmap represent some
  interval, e.g., the first block would be [floor(sqrt(n)) + 1,
  floor(sqrt(n)) + s]. Then iterate over the list of primes below sqrt(n)
  to strike all composite numbers in the block. Then output anything you
  need those primes for, and reuse the bitmap storage for the next
  interval.
  
  I think the trialdiv code currently in gmp also includes some sieve?
  
There is a prime sieve in nextprime.c at the top level.  It uses very
little memory, it doesn't even table small primes, but uses a wheel for
finding numbers with which to sieve.

The sieving area uses a fixed number of chars, currently 512.

90% of the time is spent doing integer divisions...

These divisions are used for computing the index into the sieving table,
given that index 0 represents the odd number s0.

Once the sieving number (given by the wheel) becomes greater than 512*2,
it will be used to eliminate no number half of the time, since its index
will be >= 512.

This code could be sped up be perhaps 50x by adding a remainder table of
sieving numbers, since that allows one to eliminate the divisions.

Th sieving loop (for numbers >= 512) would then consist of subtracting
2*512 from these remainders, and do a sieve write operation whenever a
number gets negative.  The negative offset is modified by adding the
sieving number.

Still, we will do lots of remainder table updates, which would
ultimately cause speed problems.  If that becomes a practical problem,
then one could sort the indexes using a minheap.

Writing this is on my TODO list.  But if you beat me to it, I am happy.

(A question is how many bits to use for representing remainders, on
64-bit machines.  24 bits should suffice.)

-- 
Torbjörn