A reference source for primes up to 64K (for unit tests)

Question

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);

Answer 1

I needed a reference source for primes up to 2¹⁶... I did not want to include the 6542 primes bodily in the source code. Apart from the bloat,

A list of 6542 primes is not necessarily "bloat". At 20 primes per line, that's 328 LOC (lines of code) with zero complexity, compared to what you ended up writing, which took only 62 LOC but has relatively high computational complexity and dependencies on the filesystem, heap allocation, and so on.

this would necessitate a test for those included primes (in case of an editing accident that trashes something)

Use source control, such as git or svn; then your test for "editing accidents" can just be "is the source code up to date". If you're worried about editing accidents in the past, take a look at the history of the source file.

and the approach would be completely impractical for bigger ranges, like primes up to 2^32.

If you want arbitrarily large lists of primes, then yes, at some point it becomes silly to hard-code them all. In that case you're going to need a reference implementation anyway. So, why bother with hard-coding the shorter list? Just use your reference implementation for both cases. This way you don't need two implementations.

---

ReferencePrimesUpTo64K's constructor should not call Console.WriteLine; that's mixing high-level business logic ("when the object can't be constructed successfully, log a message to stdout") with low-level implementation logic ("when the file can't be opened, the object can't be constructed successfully"). What you should do if you can't construct the object successfully is throw an exception; that's how you report failure in C#. Your higher-level code can then catch that exception and log a message to stdout, or try again, or switch to a different algorithm, or pop up an alert box, or terminate, or whatever high-level business logic the high-level author wants to implement.

---

But my point is that you don't need any of that filesystem logic to begin with. Computing the 6542 primes from 2 to 65521 should take you about a millisecond — quite possibly less time than it would take you to open and read your file on disk!

So let's discard all that filesystem code. That leaves us with the only interesting piece of your code: the prime-computation reference implementation.

internal static List<ushort> small_primes_up_to_64K ()
{
    const int n = ushort.MaxValue;

If you mean "65535", say "65535". Don't try to confuse the reader by making him look up what ushort.MaxValue represents. Just write 65535.

    int sqrt_n_halved = (int)(Math.Sqrt(n) - 1) >> 1, max_bit = (n - 1) >> 1;

These are constants, aren't they? If you mean "127", say "127". If you mean "32767", say "32767". Don't try to make the reader do math.

    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));

This code is much more complicated than it ought to be, and wastes O(N) memory. Have you benchmarked it, compared to the naive solution? Which, for reference, would be something like

internal static bool is_prime(int x)
{
    if (x <= 2) return (x == 2);  // (EDITED: see comments)
    if (x % 2 == 0) return false;
    const int limit = Math.Sqrt(x);
    for (int i=3; i <= limit; i += 2) {
        if (x % i == 0) return false;
    }
    return true;
}

internal static IEnumerable<int> small_primes_between(int m, int n)
{
    for (int i = m; i < n; ++i) {
        if (is_prime(i)) yield return i;
    }
}

I suspect you'll find that this simple code is approximately as fast as your complicated file-handling code — by which I mean, you'll never notice any difference in runtime in practice.

Answer 2

To get most of the stylistic issues out of the way, I ran your code through my project VSDiagnostics and it added braces, corrected identifiers according to naming conventions and changed a few things like using nameof() where applicable. I've also made some manual changes like using string interpolation instead of concatenation.

The end result was this:

using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;

public class ReferencePrimesUpTo64K
{
    public static string ReferenceFile = "reference_primes_up_to_64k.txt";
    public static ushort[] Primes;

    static ReferencePrimesUpTo64K()
    {
        if (!File.Exists(ReferenceFile))
        {
            using (var sw = new StreamWriter(ReferenceFile))
            { 
                foreach (var prime in SmallPrimesUpTo64K())
                {
                    sw.WriteLine(prime);
                }
            }

            Console.WriteLine($"# ReferencePrimesUpTo64K: file {ReferenceFile} generated.");
        }
        else if ((File.GetAttributes(ReferenceFile) & FileAttributes.ReadOnly) == 0)
        {
            Console.WriteLine($"# ReferencePrimesUpTo64K: file {ReferenceFile} not verified.");
        }

        Primes = File.ReadAllLines(ReferenceFile).Select(ushort.Parse).ToArray();
    }

    public static List<ushort> GetSmallPrimesBetween(int m, int n)
    {
        const int firstPrimePast64K = 65537;

        if (n < 2 || m > n)
        {
            return new List<ushort>();
        }

        if (n >= firstPrimePast64K)
        {
            throw new ArgumentOutOfRangeException(nameof(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> SmallPrimesUpTo64K()
    {
        const int n = ushort.MaxValue;

        int sqrtNHalved = (int)(Math.Sqrt(n) - 1) >> 1, maxBit = (n - 1) >> 1;
        var oddComposite = new bool[maxBit + 1];

        for (int i = 5 >> 1; i <= sqrtNHalved; ++i)
        {
            if (!oddComposite[i])
            {
                for (int p = (i << 1) + 1, j = p * p >> 1; j <= maxBit; j += p)
                {
                    oddComposite[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 <= maxBit; i += d, d ^= 3)
        {
            if (!oddComposite[i])
            {
                result.Add((ushort)((i << 1) + 1));
            }
        }

        return result;
    }
}

---

This is not my forte so I'm limited in the review I can offer. Something that jumps out though is this:

Primes = File.ReadAllLines(ReferenceFile).Select(ushort.Parse).ToArray();

This will be inefficient because it

1. Reads every line from the file, one by one
2. Parses every array entry, one by one

Essentially you're looping twice. I don't know how many entries you expect but if it's a sizable amount, you might want to consider using an old-fashioned for-loop instead of LINQ.

---

This might be a result of the mathematical nature of your code but try to keep variable names descriptive (mainly m & n). Then again, if it's an acceptable usage in the problem domain it's fine of course. For loops, i is often used and perfectly acceptable. p on the other hand is less common and when you mix it in with a j and some mathematical expressions, finding an issue (should there be one) is a hassle.

---

Answer 3

Use an existing, well-tested tool to generate a list:

echo 2 >primes
for i in {3..65535..2}
do factor $i
done | sed -n 's/: [^ ]*$//p' >>primes

It will take a while (almost two minutes on this old Q6600 machine), but you then have a file you can use as input, or transform into source code for your tests.

But maybe you're over-thinking your unit tests? Using a complete set of primes sounds like you're doing exhaustive tests, whereas well-written unit tests normally contain just sufficient cases to cover all the corner cases (e.g. test invalid inputs including zero and a negative number, test 1 (only 1 factor), 2 (the only even prime) and 3 (first odd prime), 4 (first true composite), 15 (first product of two distinct odd primes), 65535 (largest permitted input). Not the entire domain of 16-bit unsigned integers!