找到素数的概率(使用 miller-rabin 检验)

Probability of finding a prime (using miller-rabin test)

我已经实施了 Miller-Rabin 素数测试,并且每个函数似乎都可以单独正常工作。但是,当我尝试通过生成 70 位随机数来找到素数时,我的程序在找到通过 Miller-Rabin 测试(10 步)的数字之前平均生成超过 100000 个数字。这很奇怪,对于一个小于 70 位的随机奇数,它是质数的概率应该非常高(根据 Hadamard-de la Vallée Poussin 定理,大于 1/50)。我的代码可能有什么问题?随机数生成器是否有可能以非常低的概率抛出质数?我猜不是...非常欢迎任何帮助。

import random


def miller_rabin_rounds(n, t):
    '''Runs miller-rabin primallity test t times for n'''

    #  First find the values r and s such that 2^s * r = n - 1
    r = (n - 1) / 2
    s = 1
    while r % 2 == 0:
        s += 1
        r /= 2

    #  Run the test t times
    for i in range(t):
        a = random.randint(2, n - 1)
        y = power_remainder(a, r, n)

        if y != 1 and y != n - 1:
            #  check there is no j for which (a^r)^(2^j) = -1 (mod n)
            j = 0
            while j < s - 1 and y != n - 1:
                y = (y * y) % n
                if y == 1:
                    return False
                j += 1
            if y != n - 1:
                return False

    return True


def power_remainder(a, k, n):
    '''Computes (a^k) mod n efficiently by decomposing k into binary'''
    r = 1
    while k > 0:
        if k % 2 != 0:
            r = (r * a) % n
        a = (a * a) % n
        k //= 2
    return r


def random_odd(n):
    '''Generates a random odd number of max n bits'''
    a = random.getrandbits(n)
    if a % 2 == 0:
        a -= 1
    return a

if __name__ == '__main__':
    t = 10  # Number of Miller-Rabin tests per number
    bits = 70  # Number of bits of the random number
    a = random_odd(bits)
    count = 0
    while not miller_rabin_rounds(a, t):
        count += 1
        if count % 10000 == 0:
            print(count)
        a = random_odd(bits)
    print(a)

这在 python 2 而不是 python 3 中起作用的原因是两者处理整数除法的方式不同。在 python 2 中,3/2 = 1,而在 python 3 中,3/2=1.5.

看起来你应该在 python 3 中强制进行整数除法(而不是浮点数除法)。如果您更改代码以强制整数除法 (//),如下所示:

#  First find the values r and s such that 2^s * r = n - 1
r = (n - 1) // 2
s = 1
while r % 2 == 0:
    s += 1
    r //= 2

无论您使用什么 python 版本,您都应该看到正确的行为。