C# BigInteger Primality Class Sieve, isCarmichael, Deterministic, MillerRabin, SolovayStrassen, Fermat

BigPrimality.cs

public class BigPrimality
    {
        private const int PrimesTableLimit = 4096;
        private static readonly string CommonDir =
            Environment.GetFolderPath(Environment.SpecialFolder.MyDocuments) + "\\PrimesCache\\Primes\\";
        private readonly string                   PersistentPrimesListBigInt_FilePath = $"{CommonDir}PrimesList.bin";
        private readonly BigInteger               Three                               = new BigInteger(3);
        private readonly BigInteger               Two                                 = new BigInteger(2);
        public           int                      BitLength;
        public           int                      Candidates;
        private          RNGCryptoServiceProvider crng;
        private          LimitedList<double>      det = new LimitedList<double>(1000);
        public           double                   deta;
        public           int                      DeterministicCandidates;
        public           int                      DeterministicComposites;
        public           int                      FermatCandidates;
        public           int                      FermatComposites;
        public           int                      MillerRabinCandidates;
        public           int                      MillerRabinComposites;
        public           int                      Primes;
        private          int[]                    PrimesTable;
        public           int                      SieveComposites;
        public           int                      SolovayStrassenCandidates;
        public           int                      SolovayStrassenComposites;
        private readonly Stopwatch                sw = new Stopwatch();
        public BigPrimality()
        {
            if (File.Exists(PersistentPrimesListBigInt_FilePath))
                LoadFromFile();
            if (PrimesTable == null)
            {
                BuildPrimesTable();
                SaveToFile();
            }
        }
        public double SolovayStrassenRatio => (double) SolovayStrassenComposites / SolovayStrassenCandidates * 100.0;
        public double FermatRatio          => (double) FermatComposites / FermatCandidates * 100.0;
        public double MillerRabinRatio     => (double) MillerRabinComposites / MillerRabinCandidates * 100.0;
        public double DeterministicRatio   => (double) DeterministicComposites / DeterministicCandidates * 100.0;
        public double Ratio                => (double) Primes / Candidates * 100.0;
        public int    Iterations           { get; private set; }
        public void SaveToFile()
        {
            if (!Directory.Exists(CommonDir))
                Directory.CreateDirectory(CommonDir);
            using (var writer = new Writer(File.Open(PersistentPrimesListBigInt_FilePath, FileMode.OpenOrCreate), Encoding.Default))
            {
                writer.WriteInts(PrimesTable);
            }
        }
        public void LoadFromFile()
        {
            if (File.Exists(PersistentPrimesListBigInt_FilePath))
                using (var reader = new Reader(File.Open(PersistentPrimesListBigInt_FilePath, FileMode.Open), Encoding.Default))
                {
                    PrimesTable = reader.ReadInts();
                }
        }
        /// <summary>
        ///     Any bitlength greater than 64 is unworkable using the sieve test
        /// </summary>
        public void DoPrimalityTest(int bitlength = 32)
        {
            var    sw                  = new Stopwatch[6];
            var    sieveList           = new List<long>();
            var    millerRabinList     = new List<long>();
            var    fermatList          = new List<long>();
            var    solovayStrassenList = new List<long>();
            var    deterministicList   = new List<long>();
            var    factorsList         = new List<long>();
            var    pt                  = new List<PrimalityTest1>();
            var    pa                  = "";
            double DeterministicAVG;
            double sieveAVG;
            double millerRabinAVG;
            double fermatAVG;
            double solovayStrassenAVG;
            double FactorsAVG;
            for (var index = 0; index < sw.Length; ++index)
                sw[index] = new Stopwatch();
            var NumberOfCandidatesTried = 0;
            var confidencemr            = 1;
            var confidenceft            = 1;
            var confidencess            = 1;
            while (true)
            {
                foreach (var stopwatch in sw)
                    stopwatch.Reset();
                ++NumberOfCandidatesTried;
                var primalityTest = new PrimalityTest1();
                var mv2           = BigIntegerHelper.GetMaxValue(bitlength);
                var candidate     = NextRandomBigInt(3, mv2, bitlength);
                primalityTest.candidate = candidate;
                sw[0].Start();
                primalityTest.SievePrime = Sieve(candidate);
                sw[0].Stop();
                sw[1].Start();
                primalityTest.MillerRabinPrime = MillerRabin(candidate, confidencemr);
                sw[1].Stop();
                sw[2].Start();
                primalityTest.FermatPrime = Fermat(candidate, confidenceft);
                sw[2].Stop();
                sw[3].Start();
                primalityTest.SolovayStrassenPrime = SolovayStrassen(candidate, confidencess);
                sw[3].Stop();
                sw[4].Start();
                primalityTest.DeterministicPrime = Deterministic(candidate);
                sw[4].Stop();
                if (!primalityTest.SievePrime)
                {
                    sw[5].Start();
                    primalityTest.Factors.AddRange(GetFactors(candidate));
                    sw[5].Stop();
                }
                sieveList.Add(sw[0].ElapsedTicks);
                millerRabinList.Add(sw[1].ElapsedTicks);
                fermatList.Add(sw[2].ElapsedTicks);
                solovayStrassenList.Add(sw[3].ElapsedTicks);
                deterministicList.Add(sw[4].ElapsedTicks);
                if (!primalityTest.SievePrime)
                    factorsList.Add(sw[5].ElapsedTicks);
                sieveAVG = sieveList.Average();
                millerRabinAVG = millerRabinList.Average();
                fermatAVG = fermatList.Average();
                solovayStrassenAVG = solovayStrassenList.Average();
                DeterministicAVG = deterministicList.Average();
                if (!primalityTest.SievePrime)
                    FactorsAVG = factorsList.Average();
                if (primalityTest.SievePrime != primalityTest.MillerRabinPrime ||
                    primalityTest.SievePrime != primalityTest.FermatPrime ||
                    primalityTest.SievePrime != primalityTest.SolovayStrassenPrime)
                {
                    if (primalityTest.FermatPrime)
                        primalityTest.IsCarmichaelNumber = isCarmichael(candidate);
                    pt.Add(primalityTest);
                    pa = $"{100.0 - pt.Count / (double) NumberOfCandidatesTried * 100.0:0.00000}%";
                }
                if (primalityTest.SievePrime != primalityTest.MillerRabinPrime)
                    confidencemr++;
                if (primalityTest.SievePrime != primalityTest.FermatPrime)
                    confidenceft++;
                if (primalityTest.SievePrime != primalityTest.SolovayStrassenPrime)
                    confidencess++;
            }
        }
        private void BuildPrimesTable()
        {
            var tpt = new List<int>();
            for (var i = 3; i < PrimesTableLimit; ++i, ++i)
                if (Sieve(i))
                    tpt.Add(i);
            PrimesTable = tpt.ToArray();
        }
        /// <summary>
        ///     100% accurate. Unworkable over 64 bits.
        ///     https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
        /// </summary>
        public bool Sieve(BigInteger n)
        {
            var s = Sqrt(n);
            var a = Three;
            while (a < s)
            {
                if (n % a == 0)
                {
                    SieveComposites++;
                    return false;
                }
                a += 2;
            }
            return true;
        }
        private BigInteger Sqrt(BigInteger number)
        {
            BigInteger n = 0, p = 0;
            if (number == BigInteger.Zero)
                return BigInteger.Zero;
            var high = number >> 1;
            var low  = BigInteger.Zero;
            while (high > low + 1)
            {
                n = (high + low) >> 1;
                p = n * n;
                if (number < p)
                    high = n;
                else if (number > p)
                    low = n;
                else
                    break;
            }
            return number == p ? n : low;
        }
        /// <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 = 2)
        {
            BitLength = bitlength;
            var  candidate = new BigInteger(0);
            bool f;
            Iterations = 0;
            var mv = BigIntegerHelper.GetMaxValue(bitlength);
            var th1 = new Thread(() =>
                {
                    do
                    {
                        Iterations++;
                        candidate = NextRandomBigInt(0, mv, bitlength);
                        f = IsPrime(candidate, confidence);
                    } while (!f);
                })
                {Priority = ThreadPriority.Highest};
            th1.Start();
            th1.Join();
            return candidate;
        }
        /// <summary>
        ///     https://en.wikipedia.org/wiki/Primality_test
        ///     ~83% effective at finding composite numbers within the bit range 32 to 1024, +++?
        /// </summary>
        private bool Deterministic(BigInteger candidate)
        {
            DeterministicCandidates++;
            foreach (var p in PrimesTable)
                if (p < candidate)
                {
                    if (candidate % p == 0)
                    {
                        DeterministicComposites++;
                        return false;
                    }
                }
                else
                    break;
            return true;
        }
        public bool IsPrime(ulong bi, int confidence = 2)
        {
            return IsPrime(bi.ToBigInteger(), confidence);
        }
        public bool IsPrime(long bi, int confidence = 2)
        {
            return IsPrime(bi.ToBigInteger(), confidence);
        }
        public bool IsPrime(uint bi, int confidence = 2)
        {
            return IsPrime(bi.ToBigInteger(), confidence);
        }
        public bool IsPrime(int bi, int confidence = 2)
        {
            return IsPrime(bi.ToBigInteger(), confidence);
        }
        public void ResetCompositeCounts()
        {
            SolovayStrassenComposites = 0;
            SolovayStrassenCandidates = 0;
            FermatComposites = 0;
            FermatCandidates = 0;
            MillerRabinComposites = 0;
            MillerRabinCandidates = 0;
            DeterministicComposites = 0;
            DeterministicCandidates = 0;
            Candidates = 0;
            Primes = 0;
            BitLength = 0;
            SieveComposites = 0;
            det = new LimitedList<double>(1000);
        }
        /// <summary>
        ///     Check if a number is probably prime given a confidence level(Number of witnesses)
        /// </summary>
        public bool IsPrime(BigInteger candidate, int confidence = 2)
        {
            Candidates++;
            if (candidate == BigInteger.One)
                return false;
            if (candidate == Two)
                return true;
            if (candidate == Three)
                return true;
            if (candidate.IsEven)
                return false;
            sw.Restart();
            if (!Deterministic(candidate))
            {
                sw.Stop();
                return false;
            }
            sw.Stop();
            det.Add(sw.Elapsed.TotalMilliseconds * 1000);
            deta = det.Average();
            if (!MillerRabin(candidate, confidence))
                return false;
            Primes++;
            return true;
        }
        /// <summary>
        ///     https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
        ///     100% effective at finding composite numbers with a confidence level of two, within the bit range 32 to 1024, +++?
        /// </summary>
        private bool MillerRabin(BigInteger candidate, int confidence = 2)
        {
            MillerRabinCandidates++;
            var s = 0;
            var d = candidate - BigInteger.One;
            while ((d & 1) == 0)
            {
                d >>= 1;
                s++;
            }
            if (s == 0)
            {
                MillerRabinComposites++;
                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)
                    {
                        MillerRabinComposites++;
                        return false;
                    }
                    if (x == nmo)
                        break;
                }
                if (j == s)
                {
                    MillerRabinComposites++;
                    return false;
                }
            }
            return true;
        }
        /// <summary>
        ///     https://en.wikipedia.org/wiki/Solovay%E2%80%93Strassen_primality_test
        /// </summary>
        private bool SolovayStrassen(BigInteger candidate, int confidence = 2)
        {
            SolovayStrassenCandidates++;
            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)
                {
                    SolovayStrassenComposites++;
                    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 = 2)
        {
            FermatCandidates++;
            var pmo = candidate - 1;
            for (var i = 0; i < confidence; i++)
                if (BigInteger.ModPow(NextRandomBigInt(2, pmo), pmo, candidate) != 1)
                {
                    FermatComposites++;
                    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>
        private BigInteger NextRandomBigInt(BigInteger min, BigInteger max, int bitlength = 0)
        {
            if (crng == null)
                crng = new RNGCryptoServiceProvider();
            BigInteger n;
            var        ByteLength = 0;
            if (bitlength == 0)
                ByteLength = max.ToByteArray().Length;
            else
                ByteLength = bitlength >> 3;
            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);
            return n;
        }
    }
    public class PrimalityTest1
    {
        public BigInteger       candidate;
        public bool             DeterministicPrime;
        public List<BigInteger> Factors = new List<BigInteger>();
        public bool             FermatPrime;
        public bool             IsCarmichaelNumber;
        public bool             MillerRabinPrime;
        public bool             SievePrime;
        public bool             SolovayStrassenPrime;
        private 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)}";
        }
    }