Search code examples
primesverificationsieve-algorithm

How can one verify the proper operation of a sieve close to 2^64?


Small primes - up to about 1,000,000,000,000 - are readily available from various sources. The Prime Pages (utm.edu) have lists for the first 50 million primes, primos.mat.br goes up to 10^12, and programs like the one available at primesieve.org go even higher.

However, when it comes to numbers close to 2^64 there's only the ten primes mentioned on the page Primes just less than a power of two at primes.utm.edu and that seems to be it.

The primality test found there refuses to work on numbers that don't fit a double, others - elsewhere - fail to refuse and just print trash. The primesieve.org program refuses to work with numbers that aren't at least some 40 billion below 2^64, which doesn't exactly inspire confidence in the quality of their coding. The same result everywhere: nada, zilch, niente.

The cogs and gears of sieves start creaking around the 2^62 mark, and close to 2^64 there's hardly a cog that doesn't creak loudly threatening to break apart. Hence the need for testing the implementation is greatest where verification is most difficult, because of the scarcity/absence of reliable reference data. The primesieve.org program seems to be the only one that works at least up to 2^63 or thereabouts, but I don't trust it too much because of the above-mentioned issue.

So how then can one verify the proper operation of a sieve close to 2^64? Are there reliable lists somewhere for a million (or ten million or a hundred million) primes just below and above powers of two like 2^64, 2^63 and so on? Or are there reliable (trustworthy, verified, banged-on a lot) programs that yield such sequences or that can verify primes or lists of primes?

Once a sieve has been verified it can be used to produce handy lists with sums/checksums for loads of interesting ranges, but absent such lists the situation seems difficult...

P.S.: I determined the upper limit for the primesieve.org turbo siever to be UINT64_MAX - 10 * UINT32_MAX, or 0xFFFFFFF600000009. That means only the 10 * UINT32_MAX highest primes don't have any reference data at all so far...


Solution

  • Instead of looking for a pre-computed list, you could compare the output of your sieve to a different sieve. A good sieve, written by Tomás Oliveira e Silva, is available at http://sweet.ua.pt/tos/software/prime_sieve.html.

    Another way to test your code is by testing the primality of all numbers your sieve reports as prime (or conversely, testing the non-primality of all numbers your sieve does not report as prime). A good way to do that is the Baillie-Wagstaff test. You can find a good-quality implementation by Thomas R. Nicely at http://www.trnicely.net/misc/bpsw.html.

    You might also be interested in Jan Feitsma's tables of pseudoprimes at http://www.janfeitsma.nl/math/psp2/index, which are complete to 264.