BigIntegerPrime.cs

Generate or Test BigInteger Primes

using System;
using System.Numerics;
using System.Threading;
using System.Threading.Tasks;
public class BigIntegerPrime
{
    private readonly int _bitWidth;
    private readonly int[] _lowPrimes =
    {
        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
    };
    private readonly BigIntegerRng _rng;
    private readonly BigInteger    _twosixtyfour = "18446744073709551616".BigIntegerBase10();
    private readonly uint[]        _w0           = {2};
    private readonly uint[]        _w1           = {2, 3};
    private readonly uint[]        _w10          = {2, 3, 5, 7, 11, 13, 17, 19, 23};
    private readonly uint[]        _w11          = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    private readonly uint[]        _w12          = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41};
    private readonly uint[]        _w2           = {31, 73};
    private readonly uint[]        _w3           = {2, 3, 5};
    private readonly uint[]        _w4           = {2, 3, 5, 7};
    private readonly uint[]        _w5           = {2, 7, 61};
    private readonly uint[]        _w6           = {2, 13, 23, 1662803};
    private readonly uint[]        _w7           = {2, 3, 5, 7, 11};
    private readonly uint[]        _w8           = {2, 3, 5, 7, 11, 13};
    private readonly uint[]        _w9           = {2, 3, 5, 7, 11, 13, 17};
    private          BigInteger    _nextPrime;
    public BigIntegerPrime(int BitWidthHint)
    {
        _rng          = new BigIntegerRng(BitWidthHint);
        _rng.Unsigned = true;
        _rng.OddsOnly = true;
        _bitWidth     = BitWidthHint;
    }
    public bool IsPrime(BigInteger candidate)
    {
        if (candidate == 1)
            return false;
        if (candidate == 2 || candidate == 3 || candidate == 5)
            return true;
        if (candidate <= 1000)
            return CheckLowPrimes(candidate);
        if (!PrimeCheckM10LD(candidate))
            return false;
        if (!PrimeCheckM6(candidate))
            return false;
        if (!TrialDivision(candidate))
            return false;
        if (candidate < 2047)
        {
            if (!MillerRabin(candidate, _w0))
                return false;
        }
        else if (candidate > 2047 && candidate <= 1373653)
        {
            if (!MillerRabin(candidate, _w1))
                return false;
        }
        else if (candidate > 1373653 && candidate <= 9080191)
        {
            if (!MillerRabin(candidate, _w2))
                return false;
        }
        else if (candidate > 9080191 && candidate <= 25326001)
        {
            if (!MillerRabin(candidate, _w3))
                return false;
        }
        else if (candidate > 25326001 && candidate <= 3215031751)
        {
            if (!MillerRabin(candidate, _w4))
                return false;
        }
        else if (candidate > 3215031751 && candidate <= 4759123141)
        {
            if (!MillerRabin(candidate, _w5))
                return false;
        }
        else if (candidate > 4759123141 && candidate <= 1122004669633)
        {
            if (!MillerRabin(candidate, _w6))
                return false;
        }
        else if (candidate > 1122004669633 && candidate <= 2152302898747)
        {
            if (!MillerRabin(candidate, _w7))
                return false;
        }
        else if (candidate > 2152302898747 && candidate <= 3474749660383)
        {
            if (!MillerRabin(candidate, _w8))
                return false;
        }
        else if (candidate > 3474749660383 && candidate <= 341550071728321)
        {
            if (!MillerRabin(candidate, _w9))
                return false;
        }
        else if (candidate > 341550071728321 && candidate <= 3825123056546413051)
        {
            if (!MillerRabin(candidate, _w10))
                return false;
        }
        else if (candidate > 3825123056546413051 && candidate < _twosixtyfour)
        {
            if (!MillerRabin(candidate, _w11))
                return false;
        }
        else if (candidate > _twosixtyfour)
        {
            if (!MillerRabin(candidate, _w12))
                return false;
        }
        return true;
    }
    private bool MillerRabin(BigInteger candidate, uint[] w)
    {
        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 < w.Length; ++i)
        {
            BigInteger a;
            if (candidate > _twosixtyfour)
                a = _rng.Next(3, nmo);
            else
                a = w[i];
            var x = BigInteger.ModPow(a, d, candidate);
            if (x == 1 || x == nmo)
                continue;
            for (var r = 1; r < s; ++r)
            {
                x = BigInteger.ModPow(x, 2, candidate);
                if (x == 1)
                    return false;
                if (x == nmo)
                    break;
            }
            if (x == nmo)
                continue;
            return false;
        }
        return true;
    }
    private bool TrialDivision(BigInteger candidate)
    {
        for (var i = 0; i < _lowPrimes.Length; i++)
        {
            var p = _lowPrimes[i];
            if (i < p)
            {
                if (candidate % p != 0)
                    continue;
                return false;
            }
            break;
        }
        return true;
    }
    private static bool PrimeCheckM10LD(BigInteger n)
    {
        var d1 = (int) (n % 10);
        return d1 == 1 || d1 == 3 || d1 == 7 || d1 == 9;
    }
    private static bool PrimeCheckM6(BigInteger n)
    {
        var d1 = (int) (n % 6);
        return d1 == 1 || d1 == 5;
    }
    private bool CheckLowPrimes(BigInteger val)
    {
        foreach (var v in _lowPrimes)
            if (val == v)
                return true;
        return false;
    }
    private void NextPrime(int t)
    {
        if (t > Environment.ProcessorCount - 1)
            return;
        var MaxValue = (BigInteger.One << _bitWidth) - 1;
        while (true)
        {
            if (_nextPrime != 0)
                return;
            var n = _rng.NextFast(_bitWidth);
            if (_nextPrime != 0)
                return;
            if (n == 1)
                continue;
            if (n > MaxValue)
                continue;
            if (!IsPrime(n))
                continue;
            _nextPrime = n;
        }
    }
    public BigInteger GetPrime()
    {
        void invoke()
        {
            Parallel.Invoke(new ParallelOptions {MaxDegreeOfParallelism = Environment.ProcessorCount},
                () => NextPrime(1),
                () => NextPrime(2),
                () => NextPrime(3),
                () => NextPrime(4),
                () => NextPrime(5),
                () => NextPrime(6),
                () => NextPrime(7),
                () => NextPrime(8),
                () => NextPrime(9),
                () => NextPrime(10),
                () => NextPrime(11),
                () => NextPrime(12),
                () => NextPrime(13),
                () => NextPrime(14),
                () => NextPrime(15),
                () => NextPrime(16)
            );
        }
        _nextPrime = 0;
        var thread = new Thread(invoke) {Priority = ThreadPriority.Highest};
        thread.Start();
        thread.Join();
        return _nextPrime;
    }
}

BigIntegerRng.cs

Fast Balanced Distribution BigInteger Rng

using System;
using System.Linq;
using System.Numerics;
using System.Security.Cryptography;
[Serializable]
public class BigIntegerRng : RandomNumberGenerator
{
    private int                      _bitWidth;
    private byte[]                   _buffer;
    public  RNGCryptoServiceProvider _crng;
    private BigDecimal               _dBi;
    private BigInteger               _upperLimit;
    public BigIntegerRng(int bitWidthHint)
    {
        _crng       = new RNGCryptoServiceProvider();
        _upperLimit = (BigInteger.One << bitWidthHint) - 1;
        _dBi        = new BigDecimal(1) / _upperLimit;
        _buffer     = new byte[bitWidthHint >> 3];
        _bitWidth   = bitWidthHint;
    }
    public bool OddsOnly
    {
        get;
        set;
    }
    public bool Unsigned
    {
        get;
        set;
    }
    private BigDecimal Sample(BigDecimal weight)
    {
        BigInteger Internal()
        {
            _crng.GetBytes(_buffer);
            return new BigInteger(!Unsigned ? _buffer : _buffer.Concat(new byte[] {0}).ToArray());
        }
        var        rat  = Internal() * _dBi;
        BigDecimal tens = 10;
        if (weight != 0 && weight > .9)
        {
            if (rat < weight)
            {
                var or = rat / tens + weight;
                while (or > 1)
                {
                    tens *= 10;
                    or   =  rat / tens + weight;
                }
                return or;
            }
            return rat;
        }
        return rat;
    }
    public BigInteger Next(BigInteger minValue, BigInteger maxValue)
    {
        BestBitWidth(maxValue);
        var sa = Sample(0);
        var fi = (BigDecimal) (maxValue - minValue + minValue);
        var n  = (BigInteger) (sa * fi);
        n = !OddsOnly ? n : n | 1;
        return n < minValue ? SpecialRange(minValue, maxValue) : n;
    }
    private BigInteger SpecialRange(BigInteger minValue, BigInteger maxValue)
    {
        BigInteger res;
        var        fi     = (BigDecimal) (maxValue - minValue + minValue);
        var        weight = (BigDecimal) minValue / _upperLimit;
        do
        {
            var n = (BigInteger) (Sample(weight) * fi);
            res = !OddsOnly ? n : n | 1;
        } while (res > maxValue || res < minValue);
        return res;
    }
    private void BestBitWidth(BigInteger maxValue)
    {
        var bitWidth = maxValue.GetBitWidth();
        if (_bitWidth != bitWidth)
        {
            _upperLimit = (BigInteger.One << bitWidth) - 1;
            _dBi        = new BigDecimal(1) / _upperLimit;
            _buffer     = new byte[bitWidth >> 3];
            _bitWidth   = bitWidth;
        }
    }
    public BigInteger Next(BigInteger maxValue)
    {
        return Next(0, maxValue);
    }
    public BigInteger Next()
    {
        return Next(0, _upperLimit);
    }
    public BigDecimal NextBigDecimal()
    {
        return Sample(0) * _upperLimit;
    }
    public override void GetBytes(byte[] data)
    {
        if (data == null)
            throw new ArgumentException("The buffer cannot be null.");
        _crng.GetBytes(data);
    }
    public BigInteger NextFast(int bitWidth)
    {
        if (_bitWidth != bitWidth)
        {
            _buffer   = new byte[bitWidth >> 3];
            _bitWidth = bitWidth;
        }
        _crng.GetBytes(_buffer);
        return new BigInteger(!Unsigned ? _buffer : _buffer.Concat(new byte[] {0}).ToArray());
    }
}

FixedIntXPrimality.cs

Fixed BigInteger Primality

using System;
using System.Numerics;
using System.Threading;
public class FixedIntXPrimality
{
    private readonly FixedBigInteger _twoSixtyFourPlusIntMax = new FixedBigInteger("18446744073709551616", 64);
    private readonly uint[]          _w0                     = {2};
    private readonly uint[]          _w1                     = {2, 3};
    private readonly uint[]          _w10                    = {2, 3, 5, 7, 11, 13, 17, 19, 23};
    private readonly uint[]          _w11                    = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
    private readonly uint[]          _w12                    = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41};
    private readonly uint[]          _w2                     = {31, 73};
    private readonly uint[]          _w3                     = {2, 3, 5};
    private readonly uint[]          _w4                     = {2, 3, 5, 7};
    private readonly uint[]          _w5                     = {2, 7, 61};
    private readonly uint[]          _w6                     = {2, 13, 23, 1662803};
    private readonly uint[]          _w7                     = {2, 3, 5, 7, 11};
    private readonly uint[]          _w8                     = {2, 3, 5, 7, 11, 13};
    private readonly uint[]          _w9                     = {2, 3, 5, 7, 11, 13, 17};
    private readonly int[] LowPrimes =
    {
        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
    };
    private int                                  _bitWidth;
    private FixedBigIntegerRandomNumberGenerator _rng;
    public  int                                  CandidateTried;
    public FixedIntXPrimality(int bitWidth)
    {
        _bitWidth     = bitWidth;
        _rng          = new FixedBigIntegerRandomNumberGenerator(_bitWidth);
        _rng.OddsOnly = true;
        _rng.Unsigned = true;
    }
    public void ResetBitWidth(int bitWidth)
    {
        _bitWidth     = bitWidth;
        _rng          = new FixedBigIntegerRandomNumberGenerator(_bitWidth);
        _rng.OddsOnly = true;
        _rng.Unsigned = true;
    }
    private bool TrialDivision(FixedBigInteger candidate)
    {
        for (var i = 0; i < LowPrimes.Length; i++)
        {
            var p = LowPrimes[i];
            if (i < p)
            {
                if (candidate % p != 0)
                    continue;
                return false;
            }
            break;
        }
        return true;
    }
    private static bool PrimeCheckM10LD(FixedBigInteger n)
    {
        var d1 = (int) (n % 10);
        return d1 == 1 || d1 == 3 || d1 == 7 || d1 == 9;
    }
    private static bool PrimeCheckM6(FixedBigInteger n)
    {
        var d1 = (int) (n % 6);
        return d1 == 1 || d1 == 5;
    }
    public bool IsPrime(uint bi)
    {
        return IsPrime((FixedBigInteger) bi);
    }
    public bool IsPrime(int bi)
    {
        bi = Math.Abs(bi);
        return IsPrime((FixedBigInteger) bi);
    }
    public bool IsPrime(short bi)
    {
        bi = Math.Abs(bi);
        return IsPrime((FixedBigInteger) bi);
    }
    public bool IsPrime(ushort bi)
    {
        return IsPrime((FixedBigInteger) bi);
    }
    public bool IsPrime(byte bi)
    {
        return IsPrime((FixedBigInteger) bi);
    }
    public bool IsPrime(sbyte bi)
    {
        bi = Math.Abs(bi);
        return IsPrime((FixedBigInteger) bi);
    }
    private bool CheckLowPrimes(FixedBigInteger val)
    {
        foreach (var v in LowPrimes)
            if (val == v)
                return true;
        return false;
    }
    public bool IsPrime(FixedBigInteger candidate)
    {
        if (candidate == 1)
            return false;
        if (candidate == 2 || candidate == 3 || candidate == 5)
            return true;
        if (candidate <= 1000)
            return CheckLowPrimes(candidate);
        if (!PrimeCheckM10LD(candidate))
            return false;
        if (!PrimeCheckM6(candidate))
            return false;
        if (!TrialDivision(candidate))
            return false;
        if (candidate < 2047)
        {
            if (!MillerRabin(candidate, _w0))
                return false;
        }
        else if (candidate > 2047 && candidate <= 1373653)
        {
            if (!MillerRabin(candidate, _w1))
                return false;
        }
        else if (candidate > 1373653 && candidate <= 9080191)
        {
            if (!MillerRabin(candidate, _w2))
                return false;
        }
        else if (candidate > 9080191 && candidate <= 25326001)
        {
            if (!MillerRabin(candidate, _w3))
                return false;
        }
        else if (candidate > 25326001 && candidate <= 3215031751)
        {
            if (!MillerRabin(candidate, _w4))
                return false;
        }
        else if (candidate > 3215031751 && candidate <= 4759123141)
        {
            if (!MillerRabin(candidate, _w5))
                return false;
        }
        else if (candidate > 4759123141 && candidate <= 1122004669633)
        {
            if (!MillerRabin(candidate, _w6))
                return false;
        }
        else if (candidate > 1122004669633 && candidate <= 2152302898747)
        {
            if (!MillerRabin(candidate, _w7))
                return false;
        }
        else if (candidate > 2152302898747 && candidate <= 3474749660383)
        {
            if (!MillerRabin(candidate, _w8))
                return false;
        }
        else if (candidate > 3474749660383 && candidate <= 341550071728321)
        {
            if (!MillerRabin(candidate, _w9))
                return false;
        }
        else if (candidate > 341550071728321 && candidate <= 3825123056546413051)
        {
            if (!MillerRabin(candidate, _w10))
                return false;
        }
        else if (candidate > 3825123056546413051 && candidate < _twoSixtyFourPlusIntMax)
        {
            if (!MillerRabin(candidate, _w11))
                return false;
        }
        else if (candidate > _twoSixtyFourPlusIntMax)
        {
            if (!MillerRabin(candidate, _w12))
                return false;
        }
        return true;
    }
    private bool MillerRabin(FixedBigInteger candidate, uint[] w)
    {
        var s = 0;
        var d = candidate - FixedBigInteger.One;
        while ((d & 1) == 0)
        {
            d >>= 1;
            s++;
        }
        if (s == 0)
            return false;
        var nmo = candidate - FixedBigInteger.One;
        for (var i = 0; i < w.Length; ++i)
        {
            var x = BigInteger.ModPow(w[i], (BigInteger) d, (BigInteger) candidate);
            if (x == 1 || x == nmo)
                continue;
            for (var r = 1; r < s; ++r)
            {
                x = BigInteger.ModPow(x, 2, (BigInteger) candidate);
                if (x == 1)
                    return false;
                if (x == nmo)
                    break;
            }
            if (x == nmo)
                continue;
            return false;
        }
        return true;
    }
    public FixedBigInteger GetPrimeNumber(bool fixedWidth = false)
    {
        var bw = 0;
        if (!fixedWidth)
        {
            var MaxBytes = _bitWidth >> 3;
            var bwb      = new byte[2];
            _rng.GetBytes(bwb);
            bw = bwb[0] % MaxBytes + 1;
            if (bw < 1)
                bw = 1;
        }
        else
        {
            bw = _bitWidth >> 3;
        }
        var buffer = new byte[bw];
        var n      = new FixedBigInteger(0, 0);
        CandidateTried = 0;
        var btrd = new Thread(() =>
        {
            while (true)
            {
                CandidateTried++;
                _rng.GetBytes(buffer);
                n = new FixedBigInteger(buffer, _bitWidth) | 1;
                if (n == 1)
                    continue;
                if (n > n.MaxValue)
                    continue;
                if (IsPrime(n))
                    break;
            }
        }) {Priority = ThreadPriority.Highest};
        btrd.Start();
        btrd.Join();
        return n;
    }
}