米勒拉宾算法是一个著名的依赖于概率分析来获取正确率判断的素性检验算法，其目的是实现给定范围内的$O(\log_2 n)$的单点无记忆isprime(n)。理论上的空间复杂度算是$O(1)$的。为什么说是$O(\log_2 n)$的，主要是因为快速幂的复杂度，而二的幂由于梅塞素数分布稀疏，大概率是碰不到的，所以碰到的整个区间内大多数都是$n-1$的二的幂较小的数。所以说，这些复杂度的摊还期望就是系数较大的$O(\log_2 n)$。

都2202年了，早该普及米勒拉宾算法做单点素性检验了

算法的核心要点在米勒拉宾算法

先讲要实现什么：米勒拉宾算法的核心是选取的witness，被称为凭据。一个凭据被通过的意思是指不满足以下条件：

$$a^d\not\equiv 1\ (\text{mod}\ n)\, \text{and}\, \forall r\ (0\le r\le s-1) a^{2^rd}\not\equiv n-1\ (\text{mod}\ n) $$

其中$\ d\$是指$\ n-1=2^sd\$的除开$\ 2^s\$外的因子。

如果一个合数$\ n\$的某个证据$\ a\$通过测试，则称该证据$\ a\$为伪证。

米勒拉宾算法最核心的就是找到一系列的数，保证其不都是伪证。

当然，一定要特判$\ a\$与$\ n\$的关系，因为他们必须是互质的，否则会出现$\ 0\$。最好是直接取$\ \gcd\$是否为$\ 1\$来判断，但一般省事仅仅考虑$\ a\$是否整除$\ n\$。这个是考虑到有多个凭据基的时候做的操作，如果取凭据基数为$\ 1\$时，必须采用$\ \gcd\$来特判，否则获得$\ 0\$的结果是false也认为是true了。

以下是python的简单实现，同样的逻辑可以用C++等实现，注意到快速幂如果不用龟速乘之类的来限定到本数据类型范围内的话需要扩一倍计算长度，如int需要扩展到64位，64位需要扩展到128位。python对于pow(a,b,n)是采用快速幂的模幂，等价于$a^b\ \text{mod}\ n$。

关于$a^{2^rd}\ \text{mod}\ n$，要知道，取完$2^s$就是$n-1$，但事实上是用不到的，而且$a^{n-1}\ \text{mod}\ n$一定不是$-1$。该位置可求可不求，所以边界判断以两种方法都可，但不取的话可以省一次计算。使用$a^{n-1}\ \text{mod}\ n$来判断的算法称为费马素性检验，也就是那种过不了卡迈克尔数的算法。米勒拉宾正是对费马素性检验进行修正的算法，所以一次的计算量会比纯费马大。

def miller_rabin(n: int, w: list):
    r = n // ((n - 1) & (1 - n))
    flags = []
    for a in w:
        if n % a == 0:
            return n == a
        if pow(a, r, n) != 1:
            ard = []
            while r < n:
                ard.append(pow(a, r, n))
                r <<= 1
            flags.append(ard.count(n - 1) == 0)
        else:
            flags.append(False)
    return not any(flags)

C++版应用级模板(视情况自己选择优化)：

bool miller_rabin(int n, const initializer_list<int>& witness){
    if (n < 2) return false;
    int r = n / ((n - 1) & (1 - n));
    for (int w : witness) {
        if (n % w == 0) return n == w;
        if (mpow(w, r, n) != 1) {
            for (;r < n;r <<= 1) if (mpow(w, r, n) == n - 1) break;
            if (r >= n) return false;
        }
    }
    return true;
}

重点是经过人们长期的测试，发现一些数据范围内的最有效的一些凭据如下：

$$\begin{align} n < 1,373,653&\Rightarrow w\in\{2,3\} \\ n < 9,080,191&\Rightarrow w\in\{31,73\} \\ n < 4,759,123,141&\Rightarrow w\in\{2,7,61\} \left(\le2^{32}\right)\\ n < 2,152,302,898,747&\Rightarrow w\in\{2,3,5,7,11\} \\ n < M(>2^{64}) &\Rightarrow w\in\{2,325,9375,28178,450775,9780504,1795265022\}\left(\le2^{64}\right) \end{align} $$

在数学软件Maple中的isprime就采用$\{2,3,5,7,11\}$实现，而当数字超过后会采取另外的方法。

在以下网站中可以获取到$\ w\$的最优选取：Deterministic variants of the Miller-Rabin primality test

关于复杂度：如果只取$\{2,7,61\}$三个基的判断$2^{32}$以内的数的素性，由python代码的运算量衡量，有且仅有$O(\log_2 n)$

对于一般考察的pollard rho算法作素分解的题的范围是在$10^{13}$以内，一般选取$\{2, 1215, 34862, 574237825\}$作为凭据，可以以$O(\log_2 n)$的复杂度判断$21652684502221$以内的给定任意数的素性。

from random import randrange
from math import gcd
# 未优化, 本例子仅仅只是一种应用示例, 不能作为Pollard Rho算法的模板
def Pollard_Rho(n: int):
    if n % 4 == 0:
        return 2
    if miller_rabin(n, [2, 1215, 34862, 574237825]):
        return n
    c = randrange(1, n)
    f = lambda x:(x**2+c)%n
    t, r = f(0), f(f(0))
    while t != r:
        d = gcd(abs(t - r), n)
        if d > 1:
            return d
        t, r = f(t), f(f(r))