# Optimized ulong prime test using 6k+/-1 in parallel threads with C#

For primality testing of 64 bit ulong, I have optimized a very fast trial-by-division test using possible factors of the form 6k+/-1. For input numbers less than uint.MaxValue serial testing is performed. For really larger numbers, parallel processing is performed.

Again as optimized as it is, it is still trial-by-division and does not employ any Miller Rabin techniques.

I have 2 signature overloads: one is for long which works with all signed value type integers, and the other is for ulong which works with the unsigned value type integers, i.e. everything but BigInteger since it is not a value type.

public static bool IsPrime(long number)
{
if (number < 2) { return false; }
return IsPrime((ulong)number);
}

public static bool IsPrime(ulong number)
{
// Get the quick checks out of the way.
if (number < 2) { return false; }
// Dispense with multiples of 2 and 3.
if (number % 2 == 0) { return (number == 2); }
if (number % 3 == 0) { return (number == 3); }

// Another quick check to eliminate known composites.
// http://programmers.stackexchange.com/questions/120934/best-and-most-used-algorithm-for-finding-the-primality-of-given-positive-number/120963#120963
if (!(((number - 1) % 6 == 0) || ((number + 1) % 6 == 0)))
{
return false;
}

// Quick checks are over.  What remains is a POSSIBLE prime.
// Must iterate to determine absolute the answer.
// We loop over 1/6 of the required possible factors to check,
// but since we check twice in each iteration, we are actually
// checking 1/3 of the possible divisors.  This is an improvement
// over the typical naive test of odds only which tests 1/2
// of the factors.

// Though the whole number portion of the square root of ulong.MaxValue
// would fit in a uint, we want to cast to uint and then back to ulong.
ulong root = (ulong)(uint)Math.Sqrt(number);
// Fix Corner Case: Math.Sqrt error for really HUGE ulong.
if (root == 0) root = (ulong)uint.MaxValue;

// For small enough numbers, serial is faster than parallel.
// Obviously there is some number where parallel becomes faster.
// I do not know at which point that occurs.
// I have arbitrarily chosen a point based on the square root.
//     rootCutoff = 33,554,432
//     square of rootCutoff = 1,125,899,906,842,624
const ulong rootCutoff = 65536UL * 512UL;
if ((root < rootCutoff) || (Environment.ProcessorCount == 1))
{
// Serial Loop for smaller numbers:
// Start at 5, which is (6k-1) where k=1.
// Increment the loop by 6, which is same as incrementing k by 1.
for (ulong factor = 5; factor <= root; factor += 6)
{
// Check (6k-1)
if (number % factor == 0) { return false; }
// Check (6k+1)
if (number % (factor + 2) == 0) { return false; }
}
return true;
}

// Parallel Looping for the bigger numbers:
return IsPrimeParallel(number, root);
}

private static bool IsPrimeParallel(ulong number, ulong root)
{
int composite = 0;
// I arbitrarily choose the number of chunks to be the same as the number of processors.
// Each chunk will processed in its own thread, but this in no way is equivalent
// of saying each thread goes to its own core, or vice versa.
int chunks = Environment.ProcessorCount;
// perform cast once rather than a billion times
ulong chunks64 = (ulong)chunks;
int kEnd = (int)(root / 6) + 1;
int iEnd = (kEnd / chunks) + 1;
Parallel.For(0, chunks, (chunk, loopState) =>
{
// perform cast once rather than 90-715 million times
ulong offset = (ulong)chunk + 1;
for (int i = 0; i < iEnd; i++)
{
ulong k = (chunks64 * (ulong)i) + offset;
if (k > root) { break; }
ulong factor = (6 * k) - 1;
if (number % factor == 0) // (6k-1)
{
Interlocked.Exchange(ref composite, 1);
loopState.Stop();
}
else if (number % (factor + 2) == 0) // (6k+1)
{
Interlocked.Exchange(ref composite, 1);
loopState.Stop();
}
if (loopState.IsStopped) { break; }
}
});
return (composite == 0);
}


As I’ve seen on SO and CR, parallel processing can easily be done wrong and take much longer than simple serial. Inside my loop, there are only a few, fast calculations. It would be a performance drag to create a thread for each k (worst case is over 700 million).

This normally is a great candidate for a range Partitioner but I need more than a simple range.

Consider a simplified example where the domain of k is 1 to 100 inclusively, which is to say that kEnd is 101. I also want to generate no more than 10 threads in this example.

I do NOT want ranges like:

Range 0 is {1, 2, 3 , …, 8, 9, 10}

Range 1 is {11, 12, 13 , …, 18, 19, 20}

. . .

Range 9 is {91, 92, 93 , …, 98, 99, 100}

I want series like:

Series 0 is {1, 11, 21, …, 71, 81, 91}

Series 1 is {2, 12, 22, …, 72, 82, 92}

. . .

Series 9 is {10, 20, 30, …, 80, 90, 100}

Square Root Corner Case

There is a special corner case where Math.Sqrt(ulong.MaxValue) returns the wrong value. This is not the fault of Sqrt itself but just the nature of the (implicit) cast of ulong.MaxValue to a double since you have an integer value fully and exactly represented by 64 bits but you are squeezing them into a 64 bit floating point approximation.

This corner case is easy to detect and just as easy to correct. The key is to be aware of it in the first place.

Performance

[Edit: I erroneously listed the 31 & 32 bit times as seconds. The correct unit is milliseconds.]

On my 8-core laptop in serial only mode:

Largest 31 bit prime takes 0.17 milliseconds.

Largest 32 bit prime takes 0.23 milliseconds.

Largest 63 bit prime takes over 13 seconds.

Largest 64 bit prime takes over 18 seconds.

In parallel mode:

Largest 63 bit prime takes 6.6 seconds.

Largest 64 bit prime takes 9.28 seconds.

Questions

Being this is CR, there is always an implied question of “Do you have any constructive comments?”

While I have used Parallel.For many times before, this is my first implementation where I had to create a specifically arranged series rather than a simple range. Is this done correctly and/or could it be done better?

Other than using Miller-Rabin techniques, can this be made faster? I’ve used more threads and less threads - or really what I call chunks and assuming that each chunk gets its own thread - but on my laptop the fastest times consistently were when I used a chunk count equal to my processor count.

If I understand correctly, the (k > root) test may only trigger on a few values at the end of the loop to avoid computation. If it is indeed the case then it may be cheaper to avoid the test altogether and risk doing more work.

--

I would try to have a loop on factor and throw i and k away : no multiplications, less additions.

I quickly did a test that approximates the original loop to have an idea of the gain:

    public static ulong IsPrimeInternalLoop(ulong number, ulong root)
{
ulong ret = 0;
int composite = 0;
// I arbitrarily choose the number of chunks to be the same as the number of processors.
// Each chunk will processed in its own thread, but this in no way is equivalent
// of saying each thread goes to its own core, or vice versa.
int chunks = Environment.ProcessorCount;
// perform cast once rather than a billion times
ulong chunks64 = (ulong)chunks;
int kEnd = (int)(root / 6) + 1;
int iEnd = (kEnd / chunks) + 1;

for(int chunk = 0; chunk < chunks; ++chunk)
{
// perform cast once rather than 90-715 million times
ulong offset = (ulong)chunk + 1;
for(int i = 0; i < iEnd; i++)
{
ulong k = (chunks64 * (ulong)i) + offset;
if(k > root) { break; }
ulong factor = (6 * k) - 1;
++ret;
}
}
return ret;
}

// without i:

public static ulong IsPrimeInternalLoop2(ulong number, ulong root)
{
ulong ret = 0;
int composite = 0;
// I arbitrarily choose the number of chunks to be the same as the number of processors.
// Each chunk will processed in its own thread, but this in no way is equivalent
// of saying each thread goes to its own core, or vice versa.
int chunks = Environment.ProcessorCount;
// perform cast once rather than a billion times
ulong chunks64 = (ulong)chunks;
//int kEnd = (int)(root / 6) + 1;
ulong kEnd = (root / 6) + chunks64;

for(int chunk = 0; chunk < chunks; ++chunk)
{
for(ulong k = (ulong)chunk + 1; k <= kEnd; k += chunks64)
{
ulong factor = (6 * k) - 1;
++ret;
}
}
return ret;
}

// without k:

public static ulong IsPrimeInternalLoop3(ulong number, ulong root)
{
ulong ret = 0;
int composite = 0;
// I arbitrarily choose the number of chunks to be the same as the number of processors.
// Each chunk will processed in its own thread, but this in no way is equivalent
// of saying each thread goes to its own core, or vice versa.
int chunks = Environment.ProcessorCount;
// perform cast once rather than a billion times
ulong chunks64 = (ulong)chunks;
ulong kEnd = (root / 6) + chunks64;

for(int chunk = 0; chunk < chunks; ++chunk)
{
ulong fEnd = (6 * kEnd) - 1;
ulong fStep = (6 * chunks64);

// perform cast once rather than 90-715 million times
for(ulong factor = 6 * ((ulong)chunk + 1) - 1; factor <= fEnd; factor += fStep)
{
++ret;
}
}
return ret;
}


Note: functions 2 & 3 do a bit too many tests. It would need to be more accurate. Still, the idea is there.

I tested the above with the biggest prime on 64 bits and its neighbors:

        ulong[] values = new ulong[] { 18446744073709551557 - 1, 18446744073709551557, 18446744073709551557 + 1 };

foreach(ulong v in values)
{
ulong root = (ulong)Math.Sqrt(v);
if(root <= 0)
{
root = int.MaxValue;
}
long d0 = long.MaxValue, d1 = long.MaxValue , d2 = long.MaxValue;
ulong a = 0, b = 0, c = 0;

for(int i = 0; i < 5; i++)
{
long start = DateTime.Now.Ticks;
a = Form1.IsPrimeInternalLoop(v, root);
long stop = DateTime.Now.Ticks;
d0 = Math.Min(stop - start, d0);
}

for(int i = 0; i < 5; i++)
{
long start = DateTime.Now.Ticks;
b = Form1.IsPrimeInternalLoop2(v, root);
long stop = DateTime.Now.Ticks;
d1 = Math.Min(stop - start, d1);
}

for(int i = 0; i < 5; i++)
{
long start = DateTime.Now.Ticks;
c = Form1.IsPrimeInternalLoop3(v, root);
long stop = DateTime.Now.Ticks;
d2 = Math.Min(stop - start, d2);
}
Console.WriteLine("value: {0} root: {1}", v, root);
Console.WriteLine("times: {0} {1} {2}", d0, d1, d2);
Console.WriteLine("counts: {0} {1} {2}", a, b, c);
}


The output shows smaller times on both modifications (broadly 9%) (times) The number of factors computed shows slightly more work on functions 2 and 3 than on the original one (counts)

value: 18446744073709551556 root: 4294967296

times: 15132027 13884024 13728024

counts: 715827888 715827889 715827890

value: 18446744073709551557 root: 4294967296

times: 15288027 13884024 13728024

counts: 715827888 715827889 715827890

value: 18446744073709551558 root: 4294967296

times: 15132027 13884024 13728024

counts: 715827888 715827889 715827890

--

But this is only loop cost, and against the multithread cost, the above is nothing.

Calling loopState.IsStopped may be a bit costly compared to the computations. So you may want to reduce it's usage, for example: having an arbitrary countdown only testing for it every few iterations.

This time I ran "IsPrime" in a loop. (i7 2600K/W7 VS15 release/x64)

The original code time: ~31000000

With local countdown of 4k around IsStopped: ~25000000 (I actually underestimated the gain)

When combining the above with the loop of the function #3: ~21000000

So here is the fastest code I got. (You still need to fix the loop count)

    private static bool IsPrimeParallel(ulong number, ulong root)
{
int composite = 0;
// I arbitrarily choose the number of chunks to be the same as the number of processors.
// Each chunk will processed in its own thread, but this in no way is equivalent
// of saying each thread goes to its own core, or vice versa.
int chunks = Environment.ProcessorCount;
// perform cast once rather than a billion times

Parallel.For(0, chunks, (chunk, loopState) =>
{
// perform cast once rather than 90-715 million times
int iscd = 4096; // is stopping count down
// access of variables outside of the lambda function may be slow, depending on how the compiler does it
ulong chunks64 = (ulong)chunks;
ulong kEnd = (root / 6) + chunks64;
ulong fEnd = (6 * kEnd) - 1;
ulong fStep = (6 * chunks64);

// perform cast once rather than 90-715 million times
for(ulong factor = 6 * ((ulong)chunk + 1) - 1; factor <= fEnd; factor += fStep)
{
if(number % factor == 0) // (6k-1)
{
Interlocked.Exchange(ref composite, 1);
loopState.Stop();
}
else if(number % (factor + 2) == 0) // (6k+1)
{
Interlocked.Exchange(ref composite, 1);
loopState.Stop();
}
if(--iscd == 0)
{
if(loopState.IsStopped) { break; }
iscd = 4096;
}
}
});
return (composite == 0);
}

• Thanks for a very detailed, thought-provoking answer. It will take me some time to properly digest and evaluate. The one issue I have at this time is the concept of iscd causing an early stop. This could result in a probable prime, whereas I would prefer a truly deterministic function. – Rick Davin Jun 4 '15 at 13:08
• For timings, I use a Stopwatch. I noticed you used DateTime.Now.Ticks. A Stopwatch is strongly preferred for timings but if you must use a DateTime then I suggest using DateTime.UtcNow at the very least. – Rick Davin Jun 4 '15 at 13:12
• iscd has no effect on determinism. It only helps avoid calling a costly synchronization function too often (IsStopped). The only negative outcome is that when a thread finds a match and calls loopState.Stop(), its siblings will most likely stop a few loops (~2K in the example) later than they would without iscd. – Hurricane Jun 4 '15 at 13:24
• I just needed an approximate measure of time difference. Now or UtcNow would have no effect on the relative time ... except when going summer/winter time so : good point. And I didn't knew about StopWatch. I'll be sure to use that one next time. – Hurricane Jun 4 '15 at 13:28
• It's a local variable. It does not need atomicity. – Hurricane Jun 4 '15 at 14:04

I've credited @Hurricane with providing the correct answer. Correct is insufficient. It was excellent, original, and thought-provoking.

I did a small reworking of his answer for the critical method. Rather than using the cryptically-named iscdfor counting down, I choose a more meaningful name of stopCheckTrigger for counting up. The up versus down doesn't matter but here at CodeReview we frown upon cryptic names.

Expanding upon his notion that performance is better accessing variables local to the lamba, I added a localNumber. REMEMBER This tweaking is only for when you expect to have well in excess of half a billion loops.

For the largest 63 bit, it went from 6.5 seconds to 6.2. For the largest 64 bit prime, it went from 9.5 seconds to 5.95.

Other than that bit of reworking, its the same as @Hurricane's answer:

private static bool IsPrimeParallel(ulong number, ulong root)
{
int composite = 0;

// I arbitrarily choose the number of chunks to be the same as the number of processors.
// Each chunk will processed in its own thread, but this in no way is equivalent
// of saying each thread goes to its own core, or vice versa.
int chunks = Environment.ProcessorCount;

Parallel.For(0, chunks, (chunk, loopState) =>
{
// perform cast once rather than 90-715 million times
// access of variables outside of the lambda function may be slow, depending on how the compiler does it.
ulong chunks64 = (ulong)chunks;
ulong kEnd = (root / 6) + chunks64;
ulong fEnd = (6 * kEnd) - 1;
ulong fStep = (6 * chunks64);

// Make a copy of localNumber inside the lamba!
ulong localNumber = number;

// This declaration MUST BE local to the lamba.
int stopCheckTrigger = 0;

// perform cast once rather than 90-715 million times
for (ulong factor = 6 * ((ulong)chunk + 1) - 1; factor <= fEnd; factor += fStep)
{
if (localNumber % factor == 0) // (6k-1)
{
Interlocked.Exchange(ref composite, 1);
loopState.Stop();
break;
}
else if (localNumber % (factor + 2) == 0) // (6k+1)
{
Interlocked.Exchange(ref composite, 1);
loopState.Stop();
break;
}

const int CheckToStopTriggerPoint = 4096;

if (++stopCheckTrigger == CheckToStopTriggerPoint)
{
if (loopState.IsStopped) { break; }
stopCheckTrigger = 0; // reset to 0
}
}
});

return (composite == 0);
}


If you are a fan of primes, I strongly encourage you to vote for @Hurricane's answer.