Downloads for post:
The previous post presented a UInt32 (uint) implementation of a deterministic Miller-Rabin Primality test. This post presents a very similar algorithm, but for UInt64 (ulong) values as well as BigInteger values, although only probabilistic for the latter.
A key part of implementing the Miller-Rabin algorithm is having access to a ModPow (or PowMod) function that computes, without overflowing, the modulo value of some number raised to some power. The previous post contained an implementation using a technique that only needs to deal with intermediate values as high as the MaxValue of the type squared. That was easily achieved by using UInt64 operations for intermediate values. The Miller Rabin IsPrime function itself also used UInt64 for one of its calculations. When implementing a UInt64 ModPow and IsPrime function, there is no UInt128 type in .Net for us to use. There is however a BigInteger type in .Net we can use. Faced with the need to have BigInteger operations in both the IsPrime and ModPow functions, I decided to first just implement a standard probabilistic Miller Rabin IsPrime function as a BigInteger extension, and use the ModPow function that already exists in .Net’s BigInteger implementation. I wrote it in a way that would allow me to also use it to do deterministic primality testing of UInt64 values.
The code below is a BigInteger implementation of a Miller Rabin IsPrime function. Besides the value being tested, it takes an argument specifying the number of witness rounds to be performed (the more rounds, the more certain the primality test is, at the expense of CPU time), and also an optional set of witness values. If no witness values are provided, or the number of witness values provided are less than the specified witness count, the algorithm randomly generates the additional witness values it needs. For this, I wrote a simple BigInteger random number generator class that I won’t go into detail about here, but may cover in a future post. The real purpose of taking UInt64 witness values as optional parameters is to provide the mechanism for doing deterministic primality checks using predetermined witness values that guarantee a correct primality determination up to some maximum range. See the previous post for more detail on that concept.
// constants for common referenced BigInteger values
private static BigInteger two = 2;
private static BigInteger three = 3;
/// <summary>
/// Determines if the specified value is a prime number, deterministically for
/// values 64 bits and less, and probabilistically for values larger than 64 bits.
/// </summary>
/// <remarks>This method chooses the algorithm to use based on
/// the magnitude of the specified value. For smaller values, a
/// simple trial division algorithm is used. For larger values up to
/// and including 64 bit values, a deterministic version of the
/// Miller-Rabin algorithm is used. For values more than 64 bits,
/// a probabilistic Miller-Rabin algorithm with a default of 64
/// witness iterations is used.</remarks>
/// <param name="value">The value to be tested for primality.</param>
/// <param name="maxWitnessCount">The maximum number of witness iterations. Default is 64.</param>
/// <returns>True if the value is prime, otherwise, false.</returns>
public static bool IsPrime(this BigInteger value, int maxWitnessCount = 64) {
if (value < two) return false;
if (value <= ulong.MaxValue) return ((ulong)value).IsPrime();
return value.IsPrimeMR(maxWitnessCount);
}
/// <summary>
/// Determines if the specified value is a prime number, using the
/// probabilistic Miller-Rabin algorithm.
/// </summary>
/// <param name="value">The value to be tested for primality.</param>
/// <param name="maxWitnessCount">The maximum number of witness iterations.</param>
/// <returns>True if the value is probably prime, otherwise, false.</returns>
public static bool IsPrimeMR(this BigInteger value, int maxWitnessCount = 64) {
// take care of the simple cases of small primes and the
// common composites having those primes as factors
if (value <= BigInteger.One) return false;
if ((value % two) == BigInteger.Zero) return value == two;
if ((value % three) == BigInteger.Zero) return value == three;
if ((value % 5) == BigInteger.Zero) return value == 5;
if (((value % 7) == BigInteger.Zero) || ((value % 11) == BigInteger.Zero)
|| ((value % 13) == BigInteger.Zero) || ((value % 17) == BigInteger.Zero)
|| ((value % 19) == BigInteger.Zero) || ((value % 23) == BigInteger.Zero)
|| ((value % 29) == BigInteger.Zero) || ((value % 31) == BigInteger.Zero)
|| ((value % 37) == BigInteger.Zero) || ((value % 41) == BigInteger.Zero)
|| ((value % 43) == BigInteger.Zero)) {
return (value <= 43);
}
return InternalIsPrimeMR(value, maxWitnessCount, new ulong[0]);
}
/// <summary>
/// Determines if the specified odd, >= 3 value is a prime number, using the
/// Miller-Rabin algorithm.
/// </summary>
/// <param name="value">The value to be tested for primality.</param>
/// <param name="witnesses">The witnesses to be used.</param>
/// <returns>True if the value is prime, otherwise, false.</returns>
internal static bool InternalIsPrimeMR(BigInteger value, int witnessCount, params ulong[] witnesses) {
// compute n − 1 as (2^s)·d (where d is odd)
BigInteger valLessOne = value - BigInteger.One;
BigInteger d = valLessOne / two; // we know that value is odd and valLessOne is even, so unroll 1st iter of loop
uint s = 1;
while ((d % two) == BigInteger.Zero) {
d /= two;
s++;
}
// test value against each witness
RandomBigInteger rand = null;
for (int i = 0; i < witnessCount; i++) {
BigInteger a;
if (i < witnesses.Length) {
a = witnesses[i];
if (a >= valLessOne) {
a %= value - three;
a += three;
}
} else {
if (rand == null) rand = new RandomBigInteger(3, valLessOne);
a = rand.Next();
}
BigInteger x = BigInteger.ModPow(a, d, value);
if (x == BigInteger.One) continue;
for (uint r = 1; (r < s) && (x != valLessOne); r++) {
x = BigInteger.ModPow(x, two, value);
if (x == BigInteger.One) return false;
}
if (x != valLessOne) return false;
}
// witnesses confirm value is prime
return true;
}
The algorithm is very similar to the one for UInt32 in the previous post. First a simple check of compositeness is performed against the common low prime values. If those primes are not a factor, then the Miller Rabin algorithm is executed. Checking to see if low prime values are factors is a trade-off. When benchmarking the performance of an IsPrime function, there are two cases to consider: determining that a composite number is composite, and determining that a prime number is prime. It typically takes much longer to handle the prime input than it does to handle the composite input. Checking low prime values will improve the speed of the composite case, but will slow down the speed of the prime case. That’s because all the low prime tests will be executed, all will fail, and then the entire Miller-Rabin logic gets executed. So there’s a bit of a balancing game to play deciding how many low primes to test against in that first step.
The benefit of using this single BigInteger IsPrime implementation for both BigInteger and UInt64 prime testing is that we kill two birds with one stone. Not only do we get something that can do deterministic tests of UInt64 values, but we also can use it to do something like testing primality of a 1024 bit BigInteger value.
The code below are the UInt64 extension methods that leverage off the BigInteger function above to do deterministic checks using pre-tested witness values. Besides that, they also tie into the UInt32 implementation to check values that fall within that range. The BigInteger function is quite a bit slower than the UInt32 function. So the idea is to have the code automatically use the function that will perform the fastest for the value being tested, regardless of the type of the original value.
/// <summary>
/// Determines if the specified value is a prime number.
/// </summary>
/// <remarks>This method chooses the best algorithm to use based on
/// the magnitude of the specified value. For smaller values, a
/// simple trial division algorithm is used. For larger values, a
/// deterministic version of the Miller-Rabin algorithm is used.</remarks>
/// <param name="value">The value to be tested for primality.</param>
/// <returns>True if the value is prime, otherwise, false.</returns>
public static bool IsPrime(this ulong value) {
if (value <= uint.MaxValue) {
return ((uint)value).IsPrime();
} else {
if ((value & 1) == 0) return value == 2;
return value.IsPrimeMR();
}
}
// precomputed witnesses and the maximum value up to which they guarantee primality
private static readonly ulong[] witnesses1 = new ulong[] { 9345883071009581737 };
private const ulong witnesses1Max = 341531;
private static readonly ulong[] witnesses2 = new ulong[] { 336781006125, 9639812373923155 };
private const uint witnesses2Max = 1050535501;
private static readonly ulong[] witnesses3 = new ulong[] { 4230279247111683200, 14694767155120705706, 16641139526367750375 };
private const ulong witnesses3Max = 350269456337;
private static readonly ulong[] witnesses4 = new ulong[] { 2, 141889084524735, 1199124725622454117, 11096072698276303650 };
private const ulong witnesses4Max = 55245642489451;
private static readonly ulong[] witnesses5 = new ulong[] { 2, 4130806001517, 149795463772692060, 186635894390467037, 3967304179347715805 };
private const ulong witnesses5Max = 7999252175582851;
private static readonly ulong[] witnesses6 = new ulong[] { 2, 123635709730000, 9233062284813009, 43835965440333360, 761179012939631437, 1263739024124850375 };
private const ulong witnesses6Max = 585226005592931977;
private static readonly ulong[] witnesses7 = new ulong[] { 2, 325, 9375, 28178, 450775, 9780504, 1795265022 };
/// <summary>
/// Determines if the specified value is a prime number, using the
/// Miller-Rabin algorithm.
/// </summary>
/// <remarks>This implementation is guaranteed to deterministically
/// decide primality for all positive value inputs.</remarks>
/// <param name="value">The ulong value to be tested for primality.</param>
/// <returns>True if the value is prime, otherwise, false.</returns>
public static bool IsPrimeMR(this ulong value) {
// take care of the simple cases of small primes and the
// common composites having those primes as factors
if ((value & 1) == 0) return value == 2;
if (value == 1) return false;
if ((value % 3) == 0) return value == 3;
if ((value % 5) == 0) return value == 5;
if (((value % 7) == 0) || ((value % 11) == 0) || ((value % 13) == 0)
|| ((value % 17) == 0) || ((value % 19) == 0) || ((value % 23) == 0)
|| ((value % 29) == 0) || ((value % 31) == 0) || ((value % 37) == 0)
|| ((value % 41) == 0) || ((value % 43) == 0)) {
return (value <= 43);
}
// perform Miller-Rabin with most efficient witness array for the value being tested
BigInteger bValue = (BigInteger)value;
if (value >= witnesses6Max) return BigIntegerExtensions.InternalIsPrimeMR(bValue, 7, witnesses7);
else if (value >= witnesses5Max) return BigIntegerExtensions.InternalIsPrimeMR(bValue, 6, witnesses6);
else if (value >= witnesses4Max) return BigIntegerExtensions.InternalIsPrimeMR(bValue, 5, witnesses5);
else if (value >= witnesses3Max) return BigIntegerExtensions.InternalIsPrimeMR(bValue, 4, witnesses4);
else if (value >= witnesses2Max) return BigIntegerExtensions.InternalIsPrimeMR(bValue, 3, witnesses3);
else if (value >= witnesses1Max) return BigIntegerExtensions.InternalIsPrimeMR(bValue, 2, witnesses2);
else return BigIntegerExtensions.InternalIsPrimeMR(bValue, 1, witnesses1);
}
With this code, we can determine absolutely whether an unsigned 64 bit integer (0 to 18,446,744,073,709,551,615) is a prime number or not, with at most 7 witness iterations.Performance-wise this is not as fast as the 32 bit version in the previous post. BigInteger operations are simply more time consuming than 32 bit operations. But I’ve tested it against two other BigInteger Miller Rabin implementations available on the web, and this is faster for both composite and prime BigInteger values of various sizes. For values 64 bits and smaller, it is much faster, sometimes by as much as 100 to 1000 times when it's able to avoid BigInteger math. For BigInteger values larger than 64 bits, the differences are more modest.
