Miller-Rabin 素性檢測

參考wiki—— https://en.wikipedia.org/wiki/Miller–Rabin_primality_test

其中a的選擇：

所以base可以查表，當n不夠大時：

        private static readonly uint[][] TEST_BASE = 
        {
            new uint[] { 2 },
            new uint[] { 2, 3 },
            new uint[] { 31, 73},
            new uint[] { 2, 3, 5},
            new uint[] { 2, 3, 5, 7},
            new uint[] { 2, 7, 61},
            new uint[] { 2, 13, 23, 1662803},
            new uint[] { 2, 3, 5, 7, 11},
            new uint[] { 2, 3, 5, 7, 13},
            new uint[] { 2, 3, 5, 7, 17},
            new uint[] { 2, 3, 5, 7, 11, 13, 17, 19, 23},
            new uint[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37},
            new uint[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37},
            new uint[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41},
        };

        private static readonly double[] TEST_THREADHOLD = 
        {
            1024d, 1373653d, 9080191d, 25326001d, 3215031751d, 4759123141d,
            1122004669633d, 2152302898747d, 3474749660383d, 341550071728321d,
            3825123056546413051d, 18446744073709551616d,
            318665857834031151167461d, 3317044064679887385961981d
        };
 

將上述程式碼轉為c#

/// <summary>
/// 米勒拉賓檢測法。如果<paramref name="n"/>為素數返回true，否則返回false。
/// </summary>
public static bool MillerRabinTest(int n)
{
    if (n < 2) return false;
    if (n == 2) return true;
    if (IsEven(n)) return false;
    
    //n - 1 = 2^s * r
    var nMinusOne = (uint)(n - 1);
    var s = nMinusOne.NumberOfTrailingZeros();
    var r = nMinusOne >> s;

    var index = 0;
    while (n >= TEST_THREADHOLD[index]) ++index;

    var allbase = TEST_BASE.Skip(index).Take(1).SelectMany(array => array).ToArray();
    foreach(var @base in allbase)
    {
        var mp = MonModPow(@base, r, (uint)n);
        if (mp != 1 && mp != nMinusOne)
        {
            for(var i = 1; i <= s - 1 && mp != nMinusOne; ++i)
            {
                mp = MonModPow(mp, 2, (uint)n);
                if (mp == 1) return false;
            }

            if (mp != nMinusOne) return false;
        }
    }

    return true;
}
 

其中

private static bool IsEven(int n) => (n & 1) == 0;

public static uint MonModPow(uint value, uint exponent, uint modulo)

參考另一篇日記《【日記】Montgomery Moduluar Exponentiation》