Created
December 28, 2017 20:00
-
-
Save mjs3339/7509daecdc65531c04916ad6a40d3e44 to your computer and use it in GitHub Desktop.
C# Unsigned Long Primality Class Sieve, isCarmichael, Deterministic, MillerRabin, SolovayStrassen, Fermat
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
public static class Primality | |
{ | |
public static readonly int[] Primes4k = | |
{ | |
2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, | |
157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, | |
317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, | |
499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, | |
683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, | |
883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, | |
1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, | |
1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, | |
1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553, 1559, 1567, | |
1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, | |
1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, | |
1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, | |
2099, 2111, 2113, 2129, 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, | |
2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441, | |
2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, | |
2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, | |
2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, | |
3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209, | |
3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, | |
3391, 3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, | |
3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, | |
3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, | |
3943, 3947, 3967, 3989 | |
}; | |
public static List<PrimalityTest> ptsd = new List<PrimalityTest>(); | |
public static List<PrimalityTest> ptsm = new List<PrimalityTest>(); | |
public static List<PrimalityTest> ptsf = new List<PrimalityTest>(); | |
public static List<PrimalityTest> ptss = new List<PrimalityTest>(); | |
public static HashSet<int> primeslist = new HashSet<int>(); | |
public static Stopwatch[] sw = new Stopwatch[5]; | |
public static List<long> sieveList = new List<long>(); | |
public static List<long> deterministicList = new List<long>(); | |
public static List<long> millerRabinList = new List<long>(); | |
public static List<long> fermatList = new List<long>(); | |
public static List<long> solovayStrassenList = new List<long>(); | |
private static readonly Timer _timer = new Timer(); | |
private static string dpa = ""; | |
private static string mpa = ""; | |
private static string fpa = ""; | |
private static string spa = ""; | |
public static uint[] opq; | |
public static double sieveAVG; | |
public static double deterministicAVG; | |
public static double millerRabinAVG; | |
public static double fermatAVG; | |
public static double solovayStrassenAVG; | |
private static bool Trigger; | |
private static void _timer_Elapsed(object sender, ElapsedEventArgs e) | |
{ | |
_timer.Enabled = false; | |
Trigger = true; | |
_timer.Enabled = true; | |
} | |
public static void DoPrimalityTest() | |
{ | |
_timer.Interval = 1000.0 * 30; | |
_timer.Elapsed += _timer_Elapsed; | |
for (var index = 0; index < sw.Length; ++index) | |
sw[index] = new Stopwatch(); | |
ptsd.Clear(); | |
ptsm.Clear(); | |
ptsf.Clear(); | |
ptss.Clear(); | |
sieveList.Clear(); | |
deterministicList.Clear(); | |
millerRabinList.Clear(); | |
fermatList.Clear(); | |
solovayStrassenList.Clear(); | |
primeslist.Clear(); | |
var NumberOfCandidatesTried = 0; | |
var confidence1 = 5; //reset to 1 for initial testing | |
var confidence2 = 5; | |
var confidence3 = 5; | |
_timer.Enabled = true; | |
while (true) | |
{ | |
foreach (var stopwatch in sw) | |
stopwatch.Reset(); | |
++NumberOfCandidatesTried; | |
var primalityTest = new PrimalityTest(); | |
var candidate = Math.Abs((int)CryRng.Next()) | 1; | |
primalityTest.candidate = candidate; | |
sw[0].Start(); | |
primalityTest.SievePrime = Sieve(candidate); | |
sw[0].Stop(); | |
sw[1].Start(); | |
primalityTest.DeterministicPrime = Deterministic(candidate); | |
sw[1].Stop(); | |
sw[2].Start(); | |
primalityTest.MillerRabinPrime = MillerRabin(candidate, confidence2); | |
sw[2].Stop(); | |
sw[3].Start(); | |
primalityTest.FermatPrime = Fermat(candidate, confidence1); | |
sw[3].Stop(); | |
sw[4].Start(); | |
primalityTest.SolovayStrassenPrime = SolovayStrassen(candidate, confidence3); | |
sw[4].Stop(); | |
if (Trigger) | |
{ | |
Trigger = false; | |
_timer.Enabled = false; | |
sieveList.Add(sw[0].ElapsedTicks); | |
deterministicList.Add(sw[1].ElapsedTicks); | |
millerRabinList.Add(sw[2].ElapsedTicks); | |
fermatList.Add(sw[3].ElapsedTicks); | |
solovayStrassenList.Add(sw[4].ElapsedTicks); | |
sieveAVG = sieveList.Average(); | |
deterministicAVG = deterministicList.Average(); | |
millerRabinAVG = millerRabinList.Average(); | |
fermatAVG = fermatList.Average(); | |
solovayStrassenAVG = solovayStrassenList.Average(); | |
_timer.Enabled = true; | |
} | |
if (primalityTest.SievePrime != primalityTest.DeterministicPrime) | |
{ | |
ptsd.Add(primalityTest); | |
dpa = $"{100.0 - ptsd.Count / (double)NumberOfCandidatesTried * 100.0:0.00000}%"; | |
} | |
if (primalityTest.SievePrime != primalityTest.MillerRabinPrime) | |
{ | |
primalityTest.IsCarmichaelNumber = isCarmichael(candidate); | |
ptsm.Add(primalityTest); | |
++confidence2; | |
mpa = $"{100.0 - ptsm.Count / (double)NumberOfCandidatesTried * 100.0:0.00000}%"; | |
} | |
if (primalityTest.SievePrime != primalityTest.FermatPrime) | |
{ | |
primalityTest.IsCarmichaelNumber = isCarmichael(candidate); | |
ptsf.Add(primalityTest); | |
++confidence1; | |
fpa = $"{100.0 - ptsf.Count / (double)NumberOfCandidatesTried * 100.0:0.00000}%"; | |
} | |
if (primalityTest.SievePrime != primalityTest.SolovayStrassenPrime) | |
{ | |
primalityTest.IsCarmichaelNumber = isCarmichael(candidate); | |
ptss.Add(primalityTest); | |
++confidence3; | |
spa = $"{100.0 - ptss.Count / (double)NumberOfCandidatesTried * 100.0:0.00000}%"; | |
} | |
//if (primalityTest.SievePrime && !primeslist.Contains(candidate)) | |
//{ | |
// primeslist.Add(candidate); | |
// opq = primeslist.AsParallel().OrderBy(i => i).ToArray(); | |
//} | |
} | |
} | |
/// <summary> | |
/// 100% accurate. | |
/// https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes | |
/// </summary> | |
/// <param name="candidate">The candidate.</param> | |
/// <returns></returns> | |
private static bool Sieve(int candidate) | |
{ | |
if ((candidate & 1) == 0) | |
return false; | |
var l = (int)Math.Sqrt(candidate); | |
for (var index = 3; index <= l; index += 2) | |
if (candidate % index == 0) | |
return false; | |
return true; | |
} | |
/// <summary> | |
/// 96.3% accurate. 7 ticks/candidate | |
/// https://en.wikipedia.org/wiki/Primality_test | |
/// </summary> | |
/// <param name="candidate">The candidate.</param> | |
/// <returns></returns> | |
public static bool Deterministic(int candidate) | |
{ | |
foreach (var p in Primes4k) | |
if (p < candidate) | |
{ | |
if (candidate % p == 0) | |
return false; | |
} | |
else | |
break; | |
return true; | |
} | |
/// <summary> | |
/// 99.998% accurate at confidence=3, 14 ticks/candidate | |
/// https://en.wikipedia.org/wiki/Fermat_pseudoprime | |
/// </summary> | |
/// <param name="candidate">The candidate.</param> | |
/// <param name="confidence">The confidence.</param> | |
/// <returns></returns> | |
public static bool Fermat(int candidate, int confidence) | |
{ | |
switch (candidate) | |
{ | |
case 1: | |
return false; | |
case 2: | |
return true; | |
case 3: | |
return true; | |
} | |
var pmo = candidate - 1; | |
for (var i = 0; i < confidence; i++) | |
{ | |
var a = Math.Abs((int)CryRng.Next(2, (ulong)pmo)) | 1; | |
if (modPow(a, pmo, candidate) != 1) | |
return false; | |
} | |
return true; | |
} | |
/// <summary> | |
/// Determines whether the specified n is carmichael. | |
/// https://en.wikipedia.org/wiki/Carmichael_number | |
/// </summary> | |
/// <param name="n">The n.</param> | |
/// <returns></returns> | |
public static bool isCarmichael(this int n) | |
{ | |
var factors = false; | |
var s = Math.Sqrt(n); | |
var a = 3; | |
while (a < n) | |
{ | |
if (a > s && !factors) | |
return false; | |
if (Gcd(a, n) > 1) | |
factors = true; | |
else if (modPow(a, n - 1, n) != 1) | |
return false; | |
a += 2; | |
} | |
return true; | |
} | |
public static int GetPrime(int confidence = 2) | |
{ | |
var flag = false; | |
int ti; | |
do | |
{ | |
ti = Math.Abs((int)CryRng.Next(2, int.MaxValue)) | 1; | |
switch (ti) | |
{ | |
case 0: | |
continue; | |
case 2: | |
return ti; | |
case 3: | |
return ti; | |
default: | |
flag = IsPrime(ti, confidence = 2); | |
goto case 0; | |
} | |
} while (!flag); | |
return ti; | |
} | |
public static bool IsPrime(this int candidate,int confidence=2) | |
{ | |
if (candidate == 1) | |
return false; | |
if (candidate == 2) | |
return true; | |
if (candidate == 3) | |
return true; | |
if (candidate % 2 == 0) | |
return false; | |
if (!Deterministic(candidate)) | |
return false; | |
if (!MillerRabin(candidate, confidence)) | |
return false; | |
return true; | |
} | |
/// <summary> | |
/// 99.9999% accurate at confidence=3, 17 ticks/candidate | |
/// https://en.wikipedia.org/wiki/Solovay%E2%80%93Strassen_primality_test | |
/// </summary> | |
/// <param name="candidate">The candidate.</param> | |
/// <param name="confidence">The confidence.</param> | |
/// <returns></returns> | |
public static bool SolovayStrassen(int candidate, int confidence) | |
{ | |
var cmo = candidate - 1; | |
for (var i = 0; i < confidence; i++) | |
{ | |
var r = Math.Abs((int)CryRng.Next()) | 1; | |
var a = r % cmo + 1; | |
var jacobian = (candidate + Jacobi(a, candidate)) % candidate; | |
var mod = modPow(a, cmo >> 1, candidate); | |
if ((jacobian == 0) || (mod != jacobian)) | |
return false; | |
} | |
return true; | |
} | |
/// <summary> | |
/// https://en.wikipedia.org/wiki/Euler%E2%80%93Jacobi_pseudoprime | |
/// </summary> | |
/// <param name="a">a.</param> | |
/// <param name="b">The b.</param> | |
/// <returns></returns> | |
private static int Jacobi(int a, int b) | |
{ | |
if ((b <= 0) || (b % 2 == 0)) | |
return 0; | |
var j = 1; | |
if (a < 0) | |
{ | |
a = -a; | |
if (b % 4 == 3) | |
j = -j; | |
} | |
while (a != 0) | |
{ | |
while (a % 2 == 0) | |
{ | |
a >>= 1; | |
if ((b % 8 == 3) || (b % 8 == 5)) | |
j = -j; | |
} | |
Swap(ref a, ref b); | |
if ((a % 4 == 3) && (b % 4 == 3)) | |
j = -j; | |
a %= b; | |
} | |
return b == 1 ? j : 0; | |
} | |
private static void Swap(ref int a, ref int b) | |
{ | |
var temp = a; | |
a = b; | |
b = temp; | |
} | |
/// <summary> | |
/// 99.9999% accurate at confidence=3, 14 ticks/candidate | |
/// https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test | |
/// </summary> | |
/// <param name="candidate">The candidate.</param> | |
/// <param name="confidence">The confidence.</param> | |
/// <returns></returns> | |
private static bool MillerRabin(int candidate, int confidence) | |
{ | |
switch (candidate) | |
{ | |
case 1: | |
return false; | |
case 2: | |
return true; | |
case 3: | |
return true; | |
} | |
if ((candidate & 1) == 0) | |
return false; | |
var s = 0; | |
var d = candidate - 1; | |
while ((d & 1) == 0) | |
{ | |
d >>= 1; | |
s++; | |
} | |
if (s == 0) | |
return false; //composite | |
var nmo = candidate - 1; | |
for (var i = 0; i < confidence; ++i) | |
{ | |
var a = Math.Abs((int)CryRng.Next(2, (ulong)nmo)) | 1; | |
var x = modPow(a, d, candidate); | |
if ((x == 1) || (x == nmo)) | |
continue; | |
int j; | |
for (j = 1; j < s; ++j) | |
{ | |
x = modPow(x, 2, candidate); | |
if (x == 1) | |
return false; //composite | |
if (x == nmo) | |
break; | |
} | |
if (j == s) | |
return false; //composite | |
} | |
return true; //prime with very good probability | |
} | |
public static int modPow(int b, int e, int m) | |
{ | |
long lb = b; //convert to long to avoid overflow | |
long le = e; | |
long lm = m; | |
long result = 1; | |
while (le > 0) | |
{ | |
if ((le & 1) == 1) | |
result = result % lm * (lb % lm) % lm; | |
le >>= 1; | |
lb = lb % lm * (lb % lm) % lm; | |
} | |
return (int)result; | |
} | |
/// <summary> | |
/// Greatest Common Denominator of (a,b) | |
/// </summary> | |
/// <param name="a">a.</param> | |
/// <param name="b">The b.</param> | |
/// <returns></returns> | |
public static int Gcd(int a, int b) | |
{ | |
while (b > 0) | |
{ | |
var t = a % b; | |
a = b; | |
b = t; | |
} | |
return a; | |
} | |
public class PrimalityTest | |
{ | |
public int candidate; | |
public bool DeterministicPrime; | |
public bool FermatPrime; | |
public bool IsCarmichaelNumber; | |
public bool MillerRabinPrime; | |
public bool SievePrime; | |
public bool SolovayStrassenPrime; | |
private static string BoolToYN(bool b) | |
{ | |
return b ? "Yes" : "No"; | |
} | |
public override string ToString() | |
{ | |
return $"Candidate:{candidate}," + | |
$"Sieve:{SievePrime}," + | |
$"Deterministic:{DeterministicPrime}," + | |
$"Fermat:{FermatPrime}," + | |
$"MillerRabin:{MillerRabinPrime}," + | |
$" SolovayStrassen: {SolovayStrassenPrime}," + | |
$"IsCarmichaelNumber:{BoolToYN(IsCarmichaelNumber)}"; | |
} | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment