# Double factorial and primorial

Torbjorn Granlund tg at gmplib.org
Wed Dec 21 09:14:51 CET 2011

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

Torbjorn Granlund <tg at gmplib.org> writes:

> then one could sort the indexes using a minheap.

Clever!

That would give the following algorithm for generating primes up to n:

1. Generate a list of all primes p_k up to sqrt(n).

2. Create a minheap which for each k contains the smallest odd multiple
of p_k which is larger than sqrt(n), together with a link to p_k (can
this be optimized? I think we should be able to get away without lots
of divisions. E.g, for those k for which p_k^2 > sqrt(n), we can put
these squares into the heap rather than the *smallest* multiple.

3. Repeatedly extract and modify the smallest element of the heap. This
gives us all (odd) composite numbers in the range sqrt(n), n, and any
missing numbers are the primes, in order.

Is this faster than doing the sieving blockwise using a bitmap? We get a
very different memory access pattern. It might lose big if the heap (of
size sqrt(n)) doesn't fit in the cache. I don't know how bad locality
you typically get with a heap.

What I am suggesting is not exactly what you describe.

I am suggesting a blockwise sieving with blocksize S representing a
range of 2S oss numbers.

For primes < 2*S, we sieve as usually, i.e., good old Eratosthenes.  For
primes > 2*S we insert the squares in a minheap.

S should be small, almost surely a fraction of the L1D size.  That means
we we usually have the majority of the Pi(sqrt(n)) primes in a heap.

Representation is not clear.  Perhaps 19 bits for the heap keys, that
would alloow for the remaining 45 bits to be used by the prime.

--
Torbjörn
```