In order to code the tests for my number theory library (a collection of routines that proved handy for coding challenges) I needed a reference source for primes up to 2^16.
I did not want to reference code in other modules to get those primes because I wanted the buck to stop right there instead of passing it on, and I did not want to include the 6542 primes bodily in the source code. Apart from the bloat, this would necessitate a test for those included primes (in case of an editing accident that trashes something) and the approach would be completely impractical for bigger ranges, like primes up to 2^32.
Writing a small sieve to get the primes is no problem, but then who verifies the correctness of that sieve code? Hen and egg and all that. So I decided that the reference source would be a text file that can be verified using whatever means and write-protected. If the file doesn't exist yet then it gets generated and a message is printed to that effect, so that the external verification can take place. Also, the class keeps nagging as long as the file is not read-only, since that presumably means it hasn't been verified yet.
Note: the static function GetSmallPrimesBetween()
is the raison d'être for this class; Primes[]
is incidental. The sieve function was ripped from my coding challenge paste library; it has a few flourishes that I would not have added if I had written it from scratch for this purpose, but I left them in because there's no reason to mess with tried and trusted code.
Final note: GetSmallPrimesBetween()
takes int
parameters instead of ushort
, and it is quite lenient regarding their values, as long as they do not result in an invalid request. This is to avoid the hassle of parameter vetting and casting at the call site.
public class ReferencePrimesUpTo64K
{
public static string ReferenceFile = "reference_primes_up_to_64k.txt";
public static ushort[] Primes;
static ReferencePrimesUpTo64K ()
{
if (!System.IO.File.Exists(ReferenceFile))
{
using (var sw = new System.IO.StreamWriter(ReferenceFile))
foreach (var prime in small_primes_up_to_64K())
sw.WriteLine(prime);
Console.WriteLine("# ReferencePrimesUpTo64K: file " + ReferenceFile + " generated.");
}
else if ((System.IO.File.GetAttributes(ReferenceFile) & FileAttributes.ReadOnly) == 0)
{
Console.WriteLine("# ReferencePrimesUpTo64K: file " + ReferenceFile + " not verified.");
}
Primes = System.IO.File.ReadAllLines(ReferenceFile).Select(ushort.Parse).ToArray();
}
public static List<ushort> GetSmallPrimesBetween (int m, int n)
{
const int FIRST_PRIME_PAST_64K = 65537;
if (n < 2 || m > n)
return new List<ushort>();
if (n >= FIRST_PRIME_PAST_64K)
throw new ArgumentOutOfRangeException("n");
int tail = Array.BinarySearch(Primes, (ushort)Math.Min(Math.Max(0, m), ushort.MaxValue));
int head = Array.BinarySearch(Primes, (ushort)Math.Min(Math.Max(0, n), ushort.MaxValue));
tail ^= tail >> 31;
return Primes.Skip(tail).Take((head < 0 ? ~head : head + 1) - tail).ToList();
}
internal static List<ushort> small_primes_up_to_64K ()
{
const int n = ushort.MaxValue;
int sqrt_n_halved = (int)(Math.Sqrt(n) - 1) >> 1, max_bit = (n - 1) >> 1;
var odd_composite = new bool[max_bit + 1];
for (int i = 5 >> 1; i <= sqrt_n_halved; ++i)
if (!odd_composite[i])
for (int p = (i << 1) + 1, j = p * p >> 1; j <= max_bit; j += p)
odd_composite[j] = true;
var result = new List<ushort>() { 2, 3 }; // mod 3 stepping on top of the mod 2 wheel
for (int i = 5 >> 1, d = 1; i <= max_bit; i += d, d ^= 3)
if (!odd_composite[i])
result.Add((ushort)((i << 1) + 1));
return result;
}
}
I only started learning C# in my spare time a few weeks ago, and so I do not have the experience for vetting this code regarding things like unnecessary inefficiencies, or things that might cause problems on some platforms or with some compilers. Advice on these points would be especially welcome.
I've put a small LINQPad script up on Pastebin that vets the reference file against a file containing the first 10000 primes which can be downloaded (manually) from the prime site for all things prime - primes.utm.edu (a.k.a. The Prime Pages).
Another route for obtaining reference primes would be via gp/PARI, using the primes()
function or the forprime()
command. There's also an web interface for it, for those who don't want to download and install the program.
In the light of Quuxplusone's objection I want to emphasise again that the actual reference is the reference file. The class is merely a proxy for the file, and the sieve is only included so that the class can be used right away and the headache of finding a source of verified primes can be deferred to a more convenient time. The final verdict can only be rendered after the file has been verified, of course, at which point a failure of the sieve (unexpected though it may be) would be uncovered.
Another consideration: a reference implementation remains verified only as long as the verified binary remains unmodified - that is, until someone hits the 'build' button again, at which point there's new binary that requires verification. Add to that the fact that I've uncovered (and reported) several compiler bugs in each of the compilers I used extensively, starting with MS VC++, Borland C++ and Delphi, and a whole slew of them in FoxPro. And there are plenty more ways how things can go wrong, even without compiler bugs (starting with stale object files/libraries). In that light a local, write-protected file is a whole lot simpler and inspires a whole lot more trust.
Note: the primes up to 64K are only the lowest level a three-level operation for sieving up to 2^64, the factors of the factors as it were. And each level gets considerably more complicated to code and verify (even 'just' a reference implementation). I've already discussed a few issues regarding the verification of prime sieves over on Stack Overflow last year:
- How can one verify the proper operation of a sieve close to 2^64?
- Checksumming large swathes of prime numbers? (for verification)
Preempting an objection that hasn't been raised yet - the class could and should compare prime count, sum and product against expected values just like the file vetting script linked above does.
assert_equal("count", 6542, Primes.Length);
var sum = Primes.Select((p) => (int)p).Sum();
assert_equal("sum", 202288087, sum);
var primorial_64bit = Primes.Select((p) => (ulong)p).Aggregate((product, prime) => product * prime);
assert_equal("product", 8987519195527561682UL, primorial_64bit);