Skip to content

Instantly share code, notes, and snippets.

@mjs3339

mjs3339/BigIntegerPrimality.cs

Last active April 30, 2022 20:53
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save mjs3339/a2fb190ecfb7badc99a97a8e2b5134cd to your computer and use it in GitHub Desktop.
Save mjs3339/a2fb190ecfb7badc99a97a8e2b5134cd to your computer and use it in GitHub Desktop.
Download ZIP
C# BigInteger Primality Class
Raw
BigIntegerPrimality.cs
public class BigIntegerPrimality
{
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
};
private readonly BigInteger Ten = new BigInteger(10);
private readonly BigInteger Three = new BigInteger(3);
private RNGCryptoServiceProvider crng;
public BigInteger BigIntegerBase10FromString(string str)
{
if(str[0] == '-')
throw new Exception("Invalid numeric string. Only positive numbers are allowed.");
var number = new BigInteger();
int i;
for(i = 0; i < str.Length; i++)
if(str[i] >= '0' && str[i] <= '9')
number = number * Ten + long.Parse(str[i].ToString());
return number;
}
/// <summary>
/// 96.3% accurate. 7 ticks/candidate
/// https://en.wikipedia.org/wiki/Primality_test
/// </summary>
public bool Deterministic(BigInteger candidate)
{
foreach(var p in Primes4k)
if(p < candidate)
{
if(candidate % p == 0)
return false;
}
else
{
break;
}
return true;
}
/// <summary>
/// 100% accurate.
/// https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
/// </summary>
public bool Sieve(BigInteger candidate)
{
if(candidate < 2 || candidate > 2 && candidate.IsEven || candidate > 2 && candidate.IsPowerOfTwo)
return false;
if(candidate == 1)
return false;
if(candidate == 2)
return true;
if(candidate == 3)
return true;
if(!Deterministic(candidate))
return false;
var rtv = true;
var s = Sqrt(candidate);
Parallel.ForEach(LinqHelper.Range(Three, s, 2),
(a, state) =>
{
if(candidate % a != 0)
return;
rtv = false;
state.Stop();
});
return rtv;
}
/// <summary>
/// https://en.wikipedia.org/wiki/Newton%27s_method
/// </summary>
public BigInteger Sqrt(BigInteger n)
{
var q = BigInteger.One << ((int) BigInteger.Log(n, 2) >> 1);
var m = BigInteger.Zero;
while(BigInteger.Abs(q - m) >= 1)
{
m = q;
q = (m + n / m) >> 1;
}
return q;
}
/// <summary>
/// https://en.wikipedia.org/wiki/Carmichael_number
/// </summary>
public bool isCarmichael(BigInteger n)
{
var factors = false;
var s = Sqrt(n);
var a = Three;
var nmo = n - 1;
while(a < n)
{
if(a > s && !factors)
return false;
if(Gcd(a, n) > 1)
factors = true;
else
if(BigInteger.ModPow(a, nmo, n) != 1)
return false;
a += 2;
}
return true;
}
public List<BigInteger> GetFactors(BigInteger n)
{
var Factors = new List<BigInteger>();
var s = Sqrt(n);
var a = Three;
while(a < s)
{
if(n % a == 0)
{
Factors.Add(a);
if(a * a != n)
Factors.Add(n / a);
}
a += 2;
}
return Factors;
}
/// <summary>
/// https://en.wikipedia.org/wiki/Greatest_common_divisor
/// </summary>
private BigInteger Gcd(BigInteger a, BigInteger b)
{
while(b > BigInteger.Zero)
{
var r = a % b;
a = b;
b = r;
}
return a;
}
/// <summary>
/// https://en.wikipedia.org/wiki/Least_common_multiple
/// </summary>
private BigInteger Lcm(BigInteger a, BigInteger b)
{
return a * b / Gcd(a, b);
}
/// <summary>
/// Get a pseudo prime number of bit length with a confidence(number of witnesses) level
/// </summary>
public BigInteger GetPrime(int bitlength, int confidence = 3)
{
var candidate = BigInteger.Zero;
bool f;
var mv = BigIntegerHelper.GetMaxValue(bitlength, false);
do
{
candidate = NextRandomBigInt(0, mv, bitlength);
f = IsPrime(candidate, confidence);
} while(!f);
return candidate;
}
public BigInteger[] GetRandomPrimesArray(int arrayLength, int bitlength, int confidence = 3)
{
var ba = new BigInteger[arrayLength];
for(var i = 0; i < arrayLength; i++)
ba[i] = GetPrime(bitlength, confidence);
return ba;
}
public BigInteger[] GetProgressivePrimesArray(BigInteger start, BigInteger end, int confidence = 3)
{
var ba = new ConcurrentBag<BigInteger>();
start = start | 1;
var rg = LinqHelper.Range(start, end, 2).ToArray();
rg.AsParallel()
.WithDegreeOfParallelism((int) Math.Floor(Environment.ProcessorCount / 100.0 * 75.0))
.ForAll(a =>
{
if(MillerRabin(a, confidence))
ba.Add(a);
});
return ba.ToArray();
}
public bool IsPrime(ulong bi, int confidence = 3)
{
return IsPrime(bi.ToBigInteger(), confidence);
}
public bool IsPrime(long bi, int confidence = 3)
{
return IsPrime(bi.ToBigInteger(), confidence);
}
public bool IsPrime(uint bi, int confidence = 3)
{
return IsPrime(bi.ToBigInteger(), confidence);
}
public bool IsPrime(int bi, int confidence = 3)
{
return IsPrime(bi.ToBigInteger(), confidence);
}
/// <summary>
/// Check if a number is probably prime given a confidence level(Number of witnesses)
/// </summary>
public bool IsPrime(BigInteger candidate, int confidence = 3)
{
if(candidate < 2 || candidate > 2 && candidate.IsEven || candidate > 2 && candidate.IsPowerOfTwo)
return false;
if(candidate == 1)
return false;
if(candidate == 2)
return true;
if(candidate == 3)
return true;
if(!MillerRabin(candidate, confidence))
return false;
return true;
}
/// <summary>
/// https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
/// ~99.99% effective at finding composite numbers with a confidence level of one, within the bit range 32 to 48, +++?
/// </summary>
public bool MillerRabin(BigInteger candidate, int confidence = 3)
{
if(candidate < 2 || candidate > 2 && candidate.IsEven || candidate > 2 && candidate.IsPowerOfTwo)
return false;
if(candidate == 1)
return false;
if(candidate == 2)
return true;
if(candidate == 3)
return true;
var s = 0;
var d = candidate - BigInteger.One;
while((d & 1) == 0)
{
d >>= 1;
s++;
}
if(s == 0)
return false;
var nmo = candidate - BigInteger.One;
for(var i = 0; i < confidence; ++i)
{
var x = BigInteger.ModPow(NextRandomBigInt(2, nmo), d, candidate);
if(x == 1 || x == nmo)
continue;
int j;
for(j = 1; j < s; ++j)
{
x = BigInteger.ModPow(x, 2, candidate);
if(x == 1)
return false;
if(x == nmo)
break;
}
if(j == s)
return false;
}
return true;
}
/// <summary>
/// https://en.wikipedia.org/wiki/Solovay%E2%80%93Strassen_primality_test
/// </summary>
public bool SolovayStrassen(BigInteger candidate, int confidence = 3)
{
if(candidate < 2 || candidate > 2 && candidate.IsEven || candidate > 2 && candidate.IsPowerOfTwo)
return false;
if(candidate == 1)
return false;
if(candidate == 2)
return true;
if(candidate == 3)
return true;
var cmo = candidate - 1;
for(var i = 0; i < confidence; i++)
{
var a = NextRandomBigInt(2, cmo) % cmo + 1;
var jacobian = (candidate + Jacobi(a, candidate)) % candidate;
if(jacobian == 0 || BigInteger.ModPow(a, cmo >> 1, candidate) != jacobian)
return false;
}
return true;
}
/// <summary>
/// https://en.wikipedia.org/wiki/Euler%E2%80%93Jacobi_pseudoprime
/// </summary>
private int Jacobi(BigInteger a, BigInteger 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;
}
var temp = a;
a = b;
b = temp;
if(a % 4 == 3 && b % 4 == 3)
j = -j;
a %= b;
}
return b == 1 ? j : 0;
}
/// <summary>
/// https://en.wikipedia.org/wiki/Fermat_pseudoprime
/// Watch for carmichael numbers (false positives)
/// </summary>
public bool Fermat(BigInteger candidate, int confidence = 3)
{
if(candidate < 2 || candidate > 2 && candidate.IsEven || candidate > 2 && candidate.IsPowerOfTwo)
return false;
if(candidate == 1)
return false;
if(candidate == 2)
return true;
if(candidate == 3)
return true;
var pmo = candidate - 1;
for(var i = 0; i < confidence; i++)
if(BigInteger.ModPow(NextRandomBigInt(2, pmo), pmo, candidate) != 1)
return false;
return true;
}
/// <summary>
/// A random number that meets the minimum and maximum limits.
/// Also ensures that the number is a positive odd number.
/// </summary>
public BigInteger NextRandomBigInt(BigInteger min, BigInteger max, int bitlength = 0)
{
if(min == max)
max += 2;
if(crng == null)
crng = new RNGCryptoServiceProvider();
BigInteger n;
var ByteLength = 0;
if(bitlength == 0)
ByteLength = max.ToByteArray().Length;
else
ByteLength = bitlength >> 3;
if(ByteLength <= 0)
throw new Exception("ByteLength cannot be zero...");
var buffer = new byte[ByteLength];
var c = 0;
do
{
crng.GetBytes(buffer);
n = new BigInteger(buffer);
c++;
} while(n < min || n >= max || n.IsEven || n.Sign == -1 || n.IsPowerOfTwo);
return n;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment