Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
NumberTheory functions, including Legendre Symbol, Tonelli-Shanks, Chinese Remainder Theorem, Modular Multiplicative Inverse, Primitive Root, Multiplicative Order, Discrete Logarithm, Discrete Root, Eulers Totient Phi, Carmichael Lambda Function, LCM, and GCD.
public static class Maths
{
/// <summary>
/// Legendre Symbol returns 1 for a (nonzero) quadratic residue mod p, -1 for a non-quadratic residue (non-residue), or 0 on zero.
/// </summary>
public static int LegendreSymbol(BigInteger a, BigInteger p)
{
if (p < 2) throw new ArgumentOutOfRangeException(nameof(p), $"Parameter '{nameof(p)}' must not be < 2, but you have supplied: {p}");
if (a == 0) return 0;
if (a == 1) return 1;
int result = -1;
if (a.Mod(2) == 0)
{
result = LegendreSymbol(a / 2, p);
if (((p * p - 1) & 8) != 0) // instead of dividing by 8, shift the mask bit
{
result = -result;
}
}
else
{
result = LegendreSymbol(p.Mod(a), a);
if (((a - 1) * (p - 1) & 4) != 0) // instead of dividing by 4, shift the mask bit
{
result = -result;
}
}
return result;
}
/// <summary>
/// Finds r such that (r | m) = goal, where (r | m) is the legendre symbol
/// </summary>
public static BigInteger LegendreSymbolSearch(BigInteger start, BigInteger modulus, BigInteger goal)
{
if (goal != -1 && goal != 0 && goal != 1)
{
throw new Exception($"Parameter '{nameof(goal)}' may only be -1, 0 or 1. It was {goal}.");
}
BigInteger counter = start;
BigInteger max = counter + modulus + 1;
do
{
if (LegendreSymbol(counter, modulus) == goal)
{
return counter;
}
counter++;
}
while (counter <= max);
throw new Exception("Legendre symbol matching criteria not found.");
}
/// <summary>
/// Finds x such that x^2 ≡ n (mod p)
/// </summary>
public static BigInteger TonelliShanks(BigInteger n, BigInteger p)
{
BigInteger ls = Maths.LegendreSymbol(n, p);
if (ls != 1) // Check that n is indeed a square
{
throw new ArithmeticException($"Parameter n is not a quadratic residue, mod p. Legendre symbol = {ls}");
}
if (p.Mod(4) == 3)
{
return BigInteger.ModPow(n, ((p + 1) / 4), p);
}
// Define odd q and s such as p-1 = q * 2^s
// That is, split p-1 up into a product of two numbers:
// an odd q, and some power of 2, v, such that v = 2^s.
BigInteger q = p - 1;
BigInteger exponent = 0;
while (q.Mod(2) == 0)
{
q = q / 2;
exponent++;
}
if (exponent == 0) { throw new Exception(); }
if (exponent == 1)
{
throw new Exception("This case should have already been covered by the p mod 4 check above.");
//return BigInteger.ModPow(n, ((p+1)/4), p);
}
//Step 1. Select a non-square z such as (z | p) = -1 , and set c ≡ z^q % p
BigInteger z = Maths.LegendreSymbolSearch(0, p, -1);
BigInteger c = BigInteger.ModPow(z, q, p); // c = z^q % p
BigInteger r = BigInteger.ModPow(n, ((q + 1) / 2), p); // r = n^(q+1/2) % p
BigInteger t = BigInteger.ModPow(n, q, p); // t = n^q % p
////Step 3: loop => {
// IF t ≡ 1 OUTUPUT r, p-r .
// ELSE find, by repeated squaring, the lowest i in (0<i<m) , such as t^(2^i) ≡ 1
// LET b ≡ c^(2^(m-i-1)), r ≡ r*b, t ≡ t*b^2 , c ≡ b^2 }
BigInteger a = 0, b = 0;
BigInteger counter = 1;
BigInteger max = exponent;
while (t != 1 && counter < max) // Maths.LegendreSymbol(t,p) != 1
{
a = BigInteger.Pow(2, (int)(max - counter - 1)); // (max-counter-1)
b = BigInteger.ModPow(c, a, p);
r = BigInteger.Multiply(r, b).Mod(p);
t = BigInteger.ModPow(t * b, 2, p);
c = BigInteger.ModPow(b, 2, p);
counter++;
}
BigInteger alt = (p - r);
bool notQuadraticResidueA = Maths.LegendreSymbol(t, p) != 1;
bool notQuadraticResidueB = Maths.LegendreSymbol(t, alt) != 1;
return r; //or (p-r)
}
/// <summary>
/// Finds the modulus m such that N ≡ a[i] (mod m) for all a[i] with 0 &lt; i &lt; a.Length
/// </summary>
public static BigInteger ChineseRemainderTheorem(BigInteger[] n, BigInteger[] a)
{
BigInteger prod = n.Aggregate(BigInteger.One, (i, j) => i * j);
BigInteger p;
BigInteger sm = 0;
for (int i = 0; i < n.Length; i++)
{
p = prod / n[i];
sm += a[i] * ModularMultiplicativeInverse(p, n[i]) * p;
}
return sm % prod;
}
/// <summary>
/// Returns x such that a*x ≡ 1 (mod m)
/// </summary>
public static BigInteger ModularMultiplicativeInverse(BigInteger a, BigInteger m)
{
if (m == 1)
{
return 0;
}
BigInteger divisor;
BigInteger dividend = a;
BigInteger diff = 0;
BigInteger result = 1;
BigInteger quotient = 0;
BigInteger lastDivisor = 0;
BigInteger remainder = m;
while (dividend > 1)
{
divisor = remainder;
quotient = BigInteger.DivRem(dividend, divisor, out remainder); // Divide
dividend = divisor;
lastDivisor = diff; // The thing to divide will be the last divisor
// Update diff and result
diff = result - quotient * diff;
result = lastDivisor;
}
if (result < 0)
{
result += m; // Make result positive
}
return result;
}
/// <summary>
/// Returns x such that a*x ≡ 1 (mod m) where m is composite
/// </summary>
public static BigInteger CompositeModularInverse(BigInteger a, BigInteger m)
{
BigInteger input = a;
input = input % m;
for (BigInteger x = 1; x < m; x++)
{
if ((input * x) % m == 1)
{
return x;
}
}
return (m + 1);
}
/// <summary>
/// Returns the number of integers less than n, that are coprime to n
/// </summary>
public static BigInteger EulersTotientPhi(BigInteger n)
{
if (n < 0) { throw new ArgumentOutOfRangeException(nameof(n)); }
if (n < 3) { return 1; }
if (n == 3) { return 2; }
if (Factorization.IsProbablePrime(n)) { return n - 1; }
BigInteger totient = n;
if ((n & 1) == 0)
{
totient >>= 1;
while (((n >>= 1) & 1) == 0)
;
}
for (BigInteger i = 3; i * i <= n; i += 2)
{
if (n % i == 0)
{
totient -= totient / i;
while ((n /= i) % i == 0)
;
}
}
if (n > 1)
{
totient -= totient / n;
}
return totient;
}
/// <summary>
/// Finds a generator of the multiplicative group of integers modulo p.
/// Said differently, its finds g such that every integer a coprime to p is congruent to a power of g modulo p.
/// g^k ≡ a (mod p)
/// </summary>
public static BigInteger FindPrimitiveRoot(BigInteger p)
{
return FindPrimitiveRoot(p, 2);
}
/// <summary>
/// Finds a generator of the multiplicative group of integers modulo p.
/// Said differently, its finds g such that every integer a coprime to p is congruent to a power of g modulo p.
/// g^k ≡ a (mod p)
/// </summary>
public static BigInteger FindPrimitiveRoot(BigInteger p, int minValue)
{
// Check if p is prime or not
if (!Factorization.IsProbablePrime(p))
{
if (Factorization.GetDistinctPrimeFactors(p).Count() != 1)
{
throw new ArgumentException($"Parameter {nameof(p)} must be a prime number.", nameof(p));
}
}
BigInteger phi = Maths.EulersTotientPhi(p); // The Euler Totient function phi of a prime number p is p-1
IEnumerable<BigInteger> primeFactors = Factorization.GetDistinctPrimeFactors(phi, BigInteger.Pow(10, 7));
List<BigInteger> powersToTry = primeFactors.Select(factor => phi / factor).ToList();
for (int g = minValue; g <= phi; g++) // Check for every number from 2 to phi
{
if (powersToTry.All(n => BigInteger.ModPow(g, n, p) != 1)) // If there was no n such that g^n ≡ 1 (mod p)
{
return g; // We found our primitive root
}
}
throw new Exception($"No primitive root found for prime {p}!"); // If no primitive root found
}
/// <summary>
/// Finds n such that b^n ≡ 1 (mod m) where m is prime
/// </summary>
public static BigInteger FindMultiplicativeOrder(BigInteger b, BigInteger m)
{
// Check if p is prime or not
if (m < 2) { throw new ArgumentException($"Parameter {nameof(m)} must be greater than 1.", nameof(m)); }
if (BigInteger.GreatestCommonDivisor(b, m) != 1) { throw new ArgumentException($"Arguments {nameof(b)} and {nameof(m)} must be coprime!"); }
m = BigInteger.Abs(m);
BigInteger phi = Maths.EulersTotientPhi(m);
var primeFactorization = Factorization.GetPrimeFactorizationTuple(phi, BigInteger.Pow(10, 7));
BigInteger result = 1;
foreach (var tuple in primeFactorization)
{
BigInteger y = phi / BigInteger.Pow(tuple.Item1, (int)tuple.Item2);
int o = 0;
BigInteger n = BigInteger.ModPow(b, y, m);
while (n > 1)
{
n = BigInteger.ModPow(n, tuple.Item1, m);
o++;
}
BigInteger o1 = BigInteger.Pow(tuple.Item1, o);
o1 = o1 / BigInteger.GreatestCommonDivisor(result, o1);
result = result * o1;
}
return result;
//throw new Exception($"No multiplicative order of {m} found!"); // If no multiplicative order found
}
/// <summary>
/// Finds n such that b^n ≡ 1 (mod m) with m not necessarily prime
/// </summary>
public static BigInteger FindCompositeMultiplicativeOrder(BigInteger b, BigInteger m)
{
if (m < 2)
{
throw new ArgumentException($"Parameter {nameof(m)} must be greater than 1.", nameof(m));
}
BigInteger phi = EulersTotientPhi(m);
IEnumerable<BigInteger> primeFactors = Factorization.GetDistinctPrimeFactors(phi, BigInteger.Pow(10, 7));
List<BigInteger> powersToTry = primeFactors.Select(factor => phi / factor).ToList();
foreach (BigInteger n in powersToTry)
{
if (BigInteger.ModPow(b, n, m) == 1) // If b^n ≡ 1 (mod m)
{
return n;
}
}
throw new Exception($"No multiplicative order of {m} found!"); // If no multiplicative order found
}
/// <summary>
/// Finds x such that a^x ≡ b (mod m) where a and m are coprime
/// </summary>
public static BigInteger DiscreteLogarithm(BigInteger a, BigInteger b, BigInteger m)
{
if (a == 0)
{
return b == 0 ? 1 : -1;
}
a %= m;
b %= m;
BigInteger k = 1;
BigInteger add = 0;
BigInteger g;
while ((g = BigInteger.GreatestCommonDivisor(a, m)) > 1)
{
if (b == k)
{
return add;
}
if (b % g != 0)
{
return -1;
}
b /= g;
m /= g;
add++;
k = (k * a / g) % m;
}
BigInteger n = m.SquareRoot() + 1;
BigInteger an = 1;
for (int i = 0; i < n; ++i)
{
an = (an * a) % m;
}
BigInteger current = b;
Dictionary<BigInteger, BigInteger> vals = new Dictionary<BigInteger, BigInteger>();
for (int q = 0; q <= n; ++q)
{
if (!vals.ContainsKey(current))
{
vals.Add(current, q);
}
current = (current * a) % m;
}
current = k;
for (int p = 1; p <= n; ++p)
{
current = (current * an) % m;
if (vals.ContainsKey(current))
{
BigInteger ans = n * p - vals[current] + add;
return ans;
}
}
return -1;
}
/// <summary>
/// Finds a such that a^x ≡ b (mod m) where m is prime
/// </summary>
public static BigInteger DiscreteRoot(BigInteger x, BigInteger b, BigInteger m)
{
if (b == 0) { return 0; }
if (x == 1) { return b; }
var g = Maths.FindPrimitiveRoot(m);
var gx = Maths.PowerMod(g, x, m);
var y = Maths.DiscreteLogarithm(gx, b, m);
return Maths.PowerMod(g, y, m);
}
/// <summary>
/// Finds the smallest k such that a^k ≡ 1 (mod n) for all a coprime to n.
/// Also known as the reduced totient function because λ(n) ≤ ϕ(n).
/// It also has the following properties: If a|b then λ(a)|λ(b) and λ(lcm(a,b)) = lcm(λ(a),λ(b)).
/// </summary>
public static BigInteger CarmichaelLambdaFunction(BigInteger n)
{
var tupl = Factorization.GetPrimeFactorizationTuple(n);
var primePowers = tupl.Select(tup => BigInteger.Pow(tup.Item1, (int)tup.Item2));
if (primePowers.Count() > 1)
{
return Maths.LCM(primePowers.Select(pp => CarmichaelLambdaFunction(pp)));
}
var primePower = primePowers.First();
BigInteger result = Maths.EulersTotientPhi(primePower);
if (primePowers.Count() == 1 && primePower % 2 == 0 && primePower > 4)
{
result /= 2;
}
return result;
}
public static BigInteger LCM(IEnumerable<BigInteger> numbers)
{
return numbers.Aggregate(LCM);
}
public static BigInteger LCM(BigInteger num1, BigInteger num2)
{
BigInteger absValue1 = BigInteger.Abs(num1);
BigInteger absValue2 = BigInteger.Abs(num2);
return (absValue1 * absValue2) / BigInteger.GreatestCommonDivisor(absValue1, absValue2);
}
public static BigInteger GCD(IEnumerable<BigInteger> numbers)
{
return numbers.Aggregate(BigInteger.GreatestCommonDivisor);
}
}
public static class Factorization
{
private static BigInteger _cacheCeiling;
private static List<BigInteger> _primeCache;
private static BigInteger _cacheLargestPrimeCurrently;
private static BigInteger[] _probablePrimeCheckBases;
static Factorization()
{
_cacheCeiling = BigInteger.Pow(11, 8);
_primeCache = new List<BigInteger> { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47 };
_probablePrimeCheckBases = _primeCache.ToArray();
_cacheLargestPrimeCurrently = _primeCache.Last();
}
public static void EnsurePrimeCacheSize(BigInteger maxPrime)
{
BigInteger boundedPrimeRequest = BigInteger.Min(maxPrime, _cacheCeiling);
if (_cacheLargestPrimeCurrently < boundedPrimeRequest)
{
_primeCache = PrimeFactory.GetPrimes(boundedPrimeRequest);
_cacheLargestPrimeCurrently = boundedPrimeRequest;
}
}
public static IEnumerable<BigInteger> GetPrimeFactorization(BigInteger value, BigInteger maxValue)
{
value = BigInteger.Abs(value);
maxValue = BigInteger.Min(maxValue, value.SquareRoot() + 1);
if (value == 0) { return new BigInteger[] { 0 }; }
if (value < 10) { if (value == 0 || value == 1 || value == 2 || value == 3 || value == 5 || value == 7) { return new List<BigInteger>() { value }; } }
BigInteger primeListSize = maxValue + 5;
EnsurePrimeCacheSize(primeListSize);
if (_primeCache.Contains(value)) { return new List<BigInteger>() { value }; }
List<BigInteger> factors = new List<BigInteger>();
foreach (BigInteger prime in _primeCache)
{
while (value % prime == 0)
{
value /= prime;
factors.Add(prime);
}
if (value == 1) break;
}
if (value != 1) { factors.Add(value); }
return factors;
}
// Essentially Miller-Rabin probabilistic primality test, with caching
public static bool IsProbablePrime(BigInteger input)
{
if (input == 2 || input == 3)
{
return true;
}
if (input < 2 || input % 2 == 0)
{
return false;
}
EnsurePrimeCacheSize(input + 1);
if (input < _cacheLargestPrimeCurrently)
{
return _primeCache.Contains(input);
}
BigInteger d = input - 1;
int s = 0;
while (d % 2 == 0)
{
d >>= 1;
s += 1;
}
foreach (BigInteger a in _probablePrimeCheckBases)
{
BigInteger x = BigInteger.ModPow(a, d, input);
if (x == 1 || x == input - 1)
{
continue;
}
for (int r = 1; r < s; r++)
{
x = BigInteger.ModPow(x, 2, input);
if (x == 1)
{
return false;
}
if (x == input - 1)
{
break;
}
}
if (x != input - 1)
{
return false;
}
}
return true;
}
public static BigInteger GetNextPrime(BigInteger fromValue)
{
BigInteger result = fromValue + 1;
if (result.IsEven)
{
result += 1;
}
while (!Factorization.IsProbablePrime(result))
{
result += 2;
}
return result;
}
public static BigInteger GetPreviousPrime(BigInteger fromValue)
{
BigInteger result = fromValue.IsEven ? fromValue - 1 : fromValue - 2;
while (result > 0)
{
if (Factorization.IsProbablePrime(result))
{
return result;
}
result -= 2;
}
throw new Exception($"No primes exist between {fromValue} and zero.");
}
public static IEnumerable<BigInteger> GetDistinctPrimeFactors(BigInteger value)
{
return GetDistinctPrimeFactors(GetPrimeFactorizationTuple(value, value.SquareRoot()));
}
public static IEnumerable<BigInteger> GetDistinctPrimeFactors(BigInteger value, BigInteger maxValue)
{
return GetDistinctPrimeFactors(GetPrimeFactorizationTuple(value, maxValue));
}
public static IEnumerable<BigInteger> GetDistinctPrimeFactors(IEnumerable<Tuple<BigInteger, BigInteger>> tuple)
{
return tuple.Select(tup => tup.Item1);
}
public static IEnumerable<Tuple<BigInteger, BigInteger>> GetPrimeFactorizationTuple(BigInteger value)
{
return GetPrimeFactorizationTuple(value, value.SquareRoot());
}
public static IEnumerable<Tuple<BigInteger, BigInteger>> GetPrimeFactorizationTuple(BigInteger value, BigInteger maxValue)
{
var factorization = GetPrimeFactorization(value, maxValue);
return GetPrimeFactorizationTuple(factorization);
}
public static IEnumerable<Tuple<BigInteger, BigInteger>> GetPrimeFactorizationTuple(IEnumerable<BigInteger> factorization)
{
BigInteger lastPrime = -1;
BigInteger primeCounter = 1;
List<Tuple<BigInteger, BigInteger>> result = new List<Tuple<BigInteger, BigInteger>>();
foreach (BigInteger prime in factorization)
{
if (prime == lastPrime)
{
primeCounter += 1;
}
else if (lastPrime != -1)
{
result.Add(new Tuple<BigInteger, BigInteger>(lastPrime, primeCounter));
primeCounter = 1;
}
lastPrime = prime;
}
result.Add(new Tuple<BigInteger, BigInteger>(lastPrime, primeCounter));
if (factorization.Distinct().Count() != result.Count)
{
throw new Exception($"There is a bug in {nameof(Factorization.GetPrimeFactorizationTuple)}!");
}
return result;
}
}
public static class PrimeFactory
{
private static BigInteger _cacheLargestPrimeCurrently;
private static BigInteger _cacheCeiling;
internal static BigInteger[] _primeCache;
static PrimeFactory()
{
_cacheCeiling = BigInteger.Pow(11, 7);
_primeCache = new BigInteger[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47 };
_cacheLargestPrimeCurrently = _primeCache.Last();
}
public static void EnsurePrimeCacheSize(BigInteger maxPrime)
{
BigInteger boundedPrimeRequest = BigInteger.Min(maxPrime, _cacheCeiling);
if (_cacheLargestPrimeCurrently < boundedPrimeRequest)
{
_primeCache = PrimeFactory.GetPrimes(boundedPrimeRequest).ToArray();
_cacheLargestPrimeCurrently = _primeCache.Last();
}
}
public static bool IsPrime(BigInteger p)
{
var absP = BigInteger.Abs(p);
EnsurePrimeCacheSize(absP + 1);
return _primeCache.Contains(absP);
}
public static int GetIndexFromValue(int value)
{
if (value == -1)
{
return -1;
}
EnsurePrimeCacheSize(value + 1);
BigInteger primeValue = _primeCache.First(p => p >= value);
int index = Array.IndexOf(_primeCache, primeValue);
return index;
}
public static BigInteger GetValueFromIndex(int index)
{
while (!_primeCache.Any() || (_primeCache.Length - 1) < index)
{
EnsurePrimeCacheSize(GetApproximateNthPrime(index) + 1);
}
BigInteger value = _primeCache[index];
return value;
}
public static int GetApproximateNthPrime(int n)
{
// n*ln( n*e*ln(ln(n)) )
double approx = n * Math.Log(n * Math.E * Math.Log(Math.Log(n)));
double ceil = Math.Ceiling(approx);
return (int)ceil;
}
public static List<BigInteger> GetPrimes(BigInteger ceiling)
{
return GetPrimes(2, ceiling);
}
public static List<BigInteger> GetPrimes(BigInteger floor, BigInteger ceiling)
{
if (floor < 2)
{
floor = 2;
}
if (floor > ceiling)
{
throw new ArgumentOutOfRangeException("floor > ceiling");
}
if (ceiling < 10)
{
if (ceiling == 9 || ceiling == 8 || ceiling == 7)
{
return new List<BigInteger>() { 2, 3, 5, 7 };
}
else if (ceiling == 6 || ceiling == 5)
{
return new List<BigInteger>() { 2, 3, 5 };
}
else if (ceiling == 4 || ceiling == 3)
{
return new List<BigInteger>() { 2, 3 };
}
else
{
return new List<BigInteger>() { 2 };
}
}
if (_primeCache.Length > 0 && _primeCache.Last() >= ceiling)
{
return _primeCache.SkipWhile(i => floor > i).TakeWhile(i => ceiling >= i).ToList();
}
BigInteger counter = 0;
BigInteger counterStart = 3;
BigInteger inc;
BigInteger sqrt = 3;
BigInteger ceil = ceiling > Int32.MaxValue ? Int32.MaxValue - 2 : ceiling;
bool[] primeMembershipArray = new bool[(int)ceil + 1];
primeMembershipArray[2] = true;
// Set all odds as true
for (counter = counterStart; counter <= ceiling; counter += 2)
{
if ((counter & 1) == 1) // Check if odd
{
primeMembershipArray[(int)counter] = true;
}
}
while (sqrt * sqrt <= ceiling)
{
counter = sqrt * sqrt;
inc = sqrt + sqrt;
while (counter <= ceiling)
{
primeMembershipArray[(int)counter] = false;
counter += inc;
}
sqrt += 2;
while (!primeMembershipArray[(int)sqrt])
{
sqrt++;
}
}
var maxRange = ceiling.SquareRoot();
IEnumerable<BigInteger> primesCollection = null;
try
{
primesCollection = Maths.GetRange(2, ceiling - 1).Where(n => primeMembershipArray[(int)(n)] == true).ToList();
}
catch (Exception ex)
{
$"Exception message:\n{ex}\n\nCeiling: {ceiling}\nMax primeMembershipArray value: {primeMembershipArray.Last()}\nMax primeMembershipArray index: {primeMembershipArray.Length - 1}\n\n".Dump();
if (ex.InnerException != null)
{
$"InnerException:\n{ex.InnerException}\n\n".Dump();
}
return new List<BigInteger>();
}
if (primesCollection.Count() > _primeCache.Count())
{
_primeCache = primesCollection.ToArray();
}
List<BigInteger> result = primesCollection.Where(l => l >= floor && l <= ceiling).ToList();
return result;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment