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

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.");
}

}

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;

}

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

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:

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


I needed a reference source for primes up to 216... 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])


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.

• Your code has quite a few issues, including the fact that it has 2 as a non-prime. You may want to post it on Code Review to get the full story... And that's precisely the reason why I put tried and trusted sieve code into my class instead of writing new code. Anyway, thanks for doing such a thorough review! I'll answer more fully when I have sufficient time, Sunday at the latest (working today + tomorrow). – DarthGizka Apr 29 '16 at 11:11
• "2 as a non-prime" — Good point. That can be fixed by changing my return false to return x == 2, so I'll edit that in now. I can also totally believe that my usage of IEnumerable and yield return is incorrect and/or unidiomatic — I'm not much of a C# programmer myself — but I'm pretty confident that returning an enumerable/iterable is strictly better than returning a whole List<int>. – Quuxplusone Apr 30 '16 at 7:56
• @Quuxplusone Yes, returning IEnumerable<T> is preferred. That is because you can materialize the IEnumerable<T>, but you cannot dematerialize a List<T>. – Der Kommissar Apr 30 '16 at 8:19
• @EBrown: good point! For the small quantities involved here the benefit would be marginal, but for a source of 32-bit or 64-bit primes it would be huge. – DarthGizka Apr 30 '16 at 16:37

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.

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.");
}

}

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;

}

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])
{
}
}

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.

• Interesting, thanks! The parameter names m and n hark back to things like SPOJ's PRIME1 etc. and I found them more convenient in practice (because of their strong definition in these tasks) than names like lo and hi, start and end (where the question arises whether the end is inclusive or exclusive). first and last would be perfectly fine, of course, but I found m and n more idiomatic. – DarthGizka Apr 29 '16 at 10:16
• @DarthGizka perfectly valid -- if the domain lingo uses those terms then it's a good identifier of course. – Jeroen Vannevel Apr 29 '16 at 14:19
• Thanks again for your review - I've learnt quite a lot from it that old hands might take for granted... For example, that the code becomes a lot more readable with the right usings (if (File.Exists()) reads almost like natural language), I've learnt about nameof, and I've filed away string interpolation and the LINQ parsing expression for detailed investigation when opportunity offers, so that in future I can strike the right balance between conciseness/readability on one hand and performance on the other. Grazie mille! – DarthGizka Apr 29 '16 at 14:35

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!

• Interesting! Unfortunately my MSYS can't find factor, though. Still a nice alternative for obtaining independent reference primes, for everyone who's got a 'nixish shell handy. As regards exhaustiveness: if I can do exhaustive in a few fractions of a second then I do exhaustive, period, and keep my mental energy for dealing with the 99.9% of cases where exhausting testing is infeasible or at least not practical. In this case, obtaining all primes below 64K and handing out the interesting sub-range is a lot simpler (a lot less error-prone) than generating the primes in a sub-range directly. – DarthGizka Apr 29 '16 at 13:04
• The function GetSmallPrimesBetween is used for verifying a similar function in the code to be tested, with extensive coverage of some 'interesting' input ranges that overlap all corner cases generously. Exhaustive testing isn't possible there anymore, though (that would be more than 64K x 64K = 4G function calls). The class will also be used in generous futz tests and other scenarios, so the effort invested in creating it pays off handsomely. By the way, the code that will be tested has a lot more corner cases, e.g. the hand-off from the wheel and presieve prime table to the sieve proper. – DarthGizka Apr 29 '16 at 13:14
• For obvious reasons, a sieve cannot return any of its wheel primes or any of its presieve primes, so these need to be pulled out of thin air. Ergo, more corner cases. – DarthGizka Apr 29 '16 at 13:16
• As regards the exhaustive or corner-case testing, I think that exposes a difference in background; my usual project work is more data-structure and control algorithms rather then numerical computation, so a different approach is likely needed. With very many separate inputs, the combinatorial explosion yields thousands of tests, even with just invalid inputs, extreme values, and one or two 'typical' cases. And some of the tests take a millisecond or more, which can mean that the tests alone can add entire seconds to the build... – Toby Speight Apr 29 '16 at 14:30
• BTW, what's MSYS? Is that your shell? (or OS?) It's not something I've heard of before. – Toby Speight Apr 29 '16 at 14:31