Number Theory functions: Legendre Symbol, Eulers Criterion, Tonelli-Shanks, Chinese Remainder Theorem, Modular Multiplicative Inverse, Primitive Root, Multiplicative Order, Discrete Logarithm, Discrete Root, Eulers Totient Phi, Carmichael Lambda, Exponentiation by Squaring, Factorial, Complex Gamma Function, Nth Root, Prime Factorization, LCM, GCD.

NumberTheory.Maths.cs

using System;
using System.Linq;
using System.Numerics;
using System.Collections.Generic;

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>
        /// Returns a^((p-1)/2) (mod p)
        /// </summary>
        public static BigInteger EulersCriterion(BigInteger a, BigInteger p)
        {
            BigInteger exponent = (p - 1) >> 1; // >> right shift == /2
            BigInteger result = BigInteger.ModPow(a, exponent, p);
            return result;
        }

	/// <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;
	}
	
	/// <summary>
        /// Least common multiple of several numbers
        /// </summary>
	public static BigInteger LCM(IEnumerable<BigInteger> numbers)
	{
		return numbers.Aggregate(LCM);
	}

	/// <summary>
	/// Least common multiple of two numbers
	/// </summary>
	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);
	}
	
	/// <summary>
        /// Greatest common divisor of several numbers
        /// </summary>
	public static BigInteger GCD(IEnumerable<BigInteger> numbers)
	{
		return numbers.Aggregate(BigInteger.GreatestCommonDivisor);
	}
	
	/// <summary>
        /// Exponentiation by squaring, using arbitrarily large signed integers
        /// </summary>
        public static BigInteger Pow(BigInteger @base, BigInteger exponent)
        {
            BigInteger b = BigInteger.Abs(@base);
            BigInteger exp = BigInteger.Abs(exponent);
            BigInteger result = BigInteger.One;
            while (exp > 0)
            {
                if ((exp & 1) == 1) // If exponent is odd (&1 == %2)
                {
                    result = (result * b);
                    exp -= 1;
                    if (exp == 0) { break; }
                }

                b = (b * b);
                exp >>= 1; // exp /= 2;
            }
            return result;
        }
	
	/// <summary>
        /// Exponentiation by squaring, modulus some number (as needed)
        /// </summary>
	public static BigInteger PowerMod(BigInteger @base, BigInteger exponent, BigInteger modulus)
        {
            BigInteger result = BigInteger.One;
            while (exponent > 0)
            {
                if ((exponent & 1) == 1) // If exponent is odd
                {
                    result = (result * @base).Mod(modulus);
                    exponent -= 1;
                    if (exponent == 0) { break; }
                }

                @base = (@base * @base).Mod(modulus);
                exponent >>= 1; // exponent /= 2;
            }
            return result.Mod(modulus);
        }
	
	/// <summary>
        /// Factorial function: n! = 1 * 2 * ... * n-1 * n
        /// </summary>
        public static BigInteger Factorial(BigInteger n)
        {
            return MultiplyRange(2, n); //Enumerable.Range(2, n-1).Product();
        }

        private static BigInteger MultiplyRange(BigInteger from, BigInteger to)
        {
            var diff = to - from;
            if (diff == 1) { return from * to; }
            if (diff == 0) { return from; }

            BigInteger half = (from + to) / 2;
            return BigInteger.Multiply(MultiplyRange(from, half), MultiplyRange(half + 1, to));
        }
	
	/// <summary>
        /// Gamma function, Γ(n); An extension of the factorial function to complex arguments, defined as: Γ(n) = (n - 1)!
        /// </summary>
        public static Complex Gamma(Complex z)
        {
            if (z.Real < 0.5) // Reflection formula
            {
                return Math.PI / (Complex.Sin(Math.PI * z) * Gamma(1 - z));
            }
            else
            {
                z -= 1;
                Complex x = P[0];
                for (var i = 1; i < G + 2; i++)
                {
                    x += P[i] / (z + i);
                }
                Complex t = z + G + 0.5;
                return Complex.Sqrt(2 * Math.PI) * (Complex.Pow(t, z + 0.5)) * Complex.Exp(-t) * x;
            }
        }

        private static int G = 7;
        private static double[] P = { 0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7 };
}

NumericExtensionMethods.cs

using System.Numerics;

public static class BigIntegerExtensionMethods
{
    public static BigInteger Clone(this BigInteger number)
    {
        return new BigInteger(number.ToByteArray());
    }

    public static BigInteger Mod(this BigInteger source, BigInteger mod)
    {
        BigInteger r = (source >= mod) ? source % mod : source;
        return (r < 0) ? BigInteger.Add(r, mod) : r.Clone();
    }

    public static BigInteger SquareRoot(this BigInteger input)
    {
        if (input.IsZero)
	{
		return new BigInteger(0);
	}

        BigInteger n = new BigInteger(0);
        BigInteger p = new BigInteger(0);
        BigInteger low = new BigInteger(0);
        BigInteger high = BigInteger.Abs(input);

        while (high > low + 1)
        {
            n = (high + low) >> 1;
            p = n * n;

            if (input < p)
            {
                high = n;
            }
            else if (input > p)
            {
                low = n;
            }
            else
            {
                break;
            }
        }
        return input == p ? n : low;
    }
}

PrimeFactorization.cs

using System;
using System.Linq;
using System.Numerics;
using System.Collections.Generic;

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;
	}
}

PrimeFactory.cs

using System;
using System.Linq;
using System.Numerics;
using System.Management;
using System.Collections.Generic;

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)
        {
            PrimesPaged primesPaged = new PrimesPaged();
            return primesPaged.Skip((int)floor).Take((int)(ceiling - floor)).ToList();
        }
        
        public class PrimesPaged : IEnumerable<BigInteger>
        {
            private static readonly uint PageSize; // L1 CPU cache size in bytes
            private static readonly uint BufferBits;
            private static readonly uint BufferBitsNext;

            static PrimesPaged()
            {
                uint cacheSize = 384000;
                List<uint> cacheSizes = CPUInfo.GetCacheSizes(CPUInfo.CacheLevel.Level1);
                if (cacheSizes.Any())
                {
                    cacheSize = cacheSizes.First() * 1000;
                }

                PageSize = cacheSize; // L1 CPU cache size in bytes
                BufferBits = PageSize * 8; // in bits
                BufferBitsNext = BufferBits * 2;
            }

            private static class CPUInfo
            {
                public static List<uint> GetCacheSizes(CacheLevel level)
                {
                    ManagementClass mc = new ManagementClass("Win32_CacheMemory");
                    ManagementObjectCollection moc = mc.GetInstances();
                    List<uint> cacheSizes = new List<uint>(moc.Count);

                    cacheSizes.AddRange(moc
                      .Cast<ManagementObject>()
                      .Where(p => (ushort)(p.Properties["Level"].Value) == (ushort)level)
                      .Select(p => (uint)(p.Properties["MaxCacheSize"].Value)));

                    return cacheSizes;
                }

                public enum CacheLevel : ushort
                {
                    Level1 = 3,
                    Level2 = 4,
                    Level3 = 5,
                }
            }

            public IEnumerator<BigInteger> GetEnumerator()
            {
                return Iterator();
            }

            IEnumerator IEnumerable.GetEnumerator()
            {
                return (IEnumerator)GetEnumerator();
            }

            private static IEnumerator<BigInteger> Iterator()
            {
                IEnumerator<BigInteger> basePrimes = null;
                List<uint> basePrimesArray = new List<uint>();
                uint[] cullBuffer = new uint[PageSize / 4]; // 4 byte words

                yield return 2;

                for (var low = (BigInteger)0; ; low += BufferBits)
                {
                    for (var bottomItem = 0; ; ++bottomItem)
                    {
                        if (bottomItem < 1)
                        {
                            if (bottomItem < 0)
                            {
                                bottomItem = 0;
                                yield return 2;
                            }

                            BigInteger next = 3 + low + low + BufferBitsNext;
                            if (low <= 0)
                            { // cull very first page
                                for (int i = 0, p = 3, sqr = 9; sqr < next; i++, p += 2, sqr = p * p)
                                {
                                    if ((cullBuffer[i >> 5] & (1 << (i & 31))) == 0)
                                    {
                                        for (int j = (sqr - 3) >> 1; j < BufferBits; j += p)
                                        {
                                            cullBuffer[j >> 5] |= 1u << j;
                                        }
                                    }
                                }
                            }
                            else
                            { // cull for the rest of the pages
                                Array.Clear(cullBuffer, 0, cullBuffer.Length);

                                if (basePrimesArray.Count == 0)
                                { // inite secondar base primes stream
                                    basePrimes = Iterator();
                                    basePrimes.MoveNext();
                                    basePrimes.MoveNext();
                                    basePrimesArray.Add((uint)basePrimes.Current); // add 3 to base primes array
                                    basePrimes.MoveNext();
                                }

                                // make sure bpa contains enough base primes...
                                for (BigInteger p = basePrimesArray[basePrimesArray.Count - 1], square = p * p; square < next;)
                                {
                                    p = basePrimes.Current;
                                    basePrimes.MoveNext();
                                    square = p * p;
                                    basePrimesArray.Add((uint)p);
                                }

                                for (int i = 0, limit = basePrimesArray.Count - 1; i < limit; i++)
                                {
                                    var p = (BigInteger)basePrimesArray[i];
                                    var start = (p * p - 3) >> 1;

                                    // adjust start index based on page lower limit...
                                    if (start >= low)
                                    {
                                        start -= low;
                                    }
                                    else
                                    {
                                        var r = (low - start) % p;
                                        start = (r != 0) ? p - r : 0;
                                    }
                                    for (var j = (uint)start; j < BufferBits; j += (uint)p)
                                    {
                                        cullBuffer[j >> 5] |= 1u << ((int)j);
                                    }
                                }
                            }
                        }
                        while (bottomItem < BufferBits && (cullBuffer[bottomItem >> 5] & (1 << (bottomItem & 31))) != 0)
                        {
                            ++bottomItem;
                        }
                        if (bottomItem < BufferBits)
                        {
                            yield return 3 + (((BigInteger)bottomItem + low) << 1);
                        }
                        else break; // outer loop for next page segment...
                    }
                }
            }
        }	
}

How do you implement BigInteger SquareRoot() function?

Edit:
Found a method in BigIntegerExtension.cs file to calculate the nth root.

But I am not able to find a function Maths.GetRange()

@gunpal5
Thanks for bringing this to my attention; there were a couple missing methods.
Here is what I changed:

• Added file: NumericExtensionMethods.cs containing the two missing methods: SquareRoot and Mod.
• Added Method: Maths.PowerMod
• Changed the PrimeFactory.GetPrimes implementation to be much more performant. This eliminates the need for Maths.GetRange. If you still need a range for whatever reason, use Enumerable.Range.