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);
}
Avoid Too Many Threads
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.
