C# 判定素数 高效算法 拉宾米勒算法

C#版本 拉宾米勒算法，这是目前最快的算法，没有之一，就是最快，9位的素数，也是秒出结果。

     private bool RabbinMillerTest(long n)
        {
            if (n < 2)
            { // 小于2的数即不是合数也不是素数
                return false;
            }
            int[] aPrimeList = {
               2, 3, 5, 7, 11, 17, 19, 23, 29, 31, 41,
                   43, 47, 53, 59, 67, 71, 73, 79, 83, 89, 97 };
            int nPrimeListSize = aPrimeList.Count();//求素数表元素个数
            for (int i = 0; i < nPrimeListSize; i++)
            {// 按照素数表中的数对当前素数进行判断
                if (n / 2 + 1 <= aPrimeList[i])
                {// 如果已经小于当前素数表的数，则一定是素数
                    return true;
                }
                if (0 == n % aPrimeList[i])
                {// 余数为0则说明一定不是素数
                    return false;
                }
            }
            // 找到r和m，使得n = 2^r * m + 1;
            long r = 0, m = n - 1; // ( n - 1 ) 一定是合数
            while (0 == (m & 1))
            {
                m >>= 1; // 右移一位
                r++; // 统计右移的次数
            }
            int nTestCnt = 10000; // 表示进行测试的次数
            Random randkey = new Random(57);
            for (int i = 0; i < nTestCnt; i++)
            { // 利用随机数进行测试，
                int a = aPrimeList[randkey.Next(0, nPrimeListSize)];
                if (1 != Montgomery(a, m, n))
                {
                    int j = 0;
                    int e = 1;
                    for (; j < r; j++)
                    {
                        if (n - 1 == Montgomery(a, m * e, n))
                        {
                            break;
                        }
                        e <<= 1;
                    }
                    if (j == r)
                    {
                        return false;
                    }
                }
            }
            return true;
        }

        long Montgomery(long n, long  p, long  m)
        { //快速计算(n^e)%m的值
            long k = 1;
            n %= m;
            while (p != 1)
            {
                if (0 != (p & 1)) k = (k * n) % m;
                n = (n * n) % m;
                p >>= 1;
            }
            return (n * k) % m;
        }