素数相关算法 Python 实现

目录

1. 定义&定理
2. 算法
   1. 素数判定
      1. 定义法
      2. Rabin-Miller 测试
   2. 素数生成
      1. 素数筛法（生成小素数列表）
      2. 生成大素数
   3. 欧几里得算法（求两数最大公因子）
   4. 扩展欧几里得算法

[TOC]

定义&定理

素数的定义

如果一个整数只有 1 和它本身两个正因子，这个数就是素数。
例如 2，3，5，7，9，13。

因子的定义

如果 b 除以 a 余数为 0，那么 a 是 b 的因子。

定理一

有 a、b、c 三个整数，如果 a 是 b 的因子，b 是 c 的因子，那么 a 是 c 的因子

定理二

有整数 n，如果 p 是 n 除 1 外的最小正因子，那么 p 是素数。

证明

​ 假设 p 不是素数，那么 p 必有一个除 1 和 p 以外的因子，记为 a，根据定理一，a 也是 n 的因子，而且 a 小于 p，与题设矛盾。

定理三

素数有无数个。

证明

​ 假设素数的个数是有限的，那么存在一个列表
$$
p_1,p_2,p_3 … p_n
$$
​ 包含所有的素数，这里的 n 就是素数的个数，那么定义
$$
P = p_1 \times p_2 \times p_3 \times … \times p_n + 1
$$
​ 取 P 除 1 外最小正因子，记为 p，根据定理二，p 是素数。但是 p 不再刚才定义的素数列表中，因为如果 p 在列表中，P 除以 p 无法整除，而会余 1。产生矛盾。

算法

素数判定

定义法

import math
def is_prime(n):
    for i in range(2, int(math.sqrt(n))+1):
        if n % i == 0:
            return False
    return True

Rabin-Miller 测试

import random

def fast_pow(a, b, mod):
    if b == 0: return 1
    if b == 1: return a % mod
    if b == 2: return a * a % mod
    if b % 2 == 0:
        tmp = fast_pow(a, b >> 1, mod)
        return tmp * tmp % mod
    return fast_pow(a, b - 1, mod) * a % mod
    
def rabin_miller(n):
    s, t = n - 1, 0
    while s % 2 == 0:
        s, t = s >> 1, t + 1
    for i in range(0, 128, 2):
        a = random.randint(2, n - 1)
        # v = a ** s % n
        v = fast_pow(a, s, n)
        print(i)
        if v != 1:
            i = 0
            while v != n - 1:
                if i == t-1:
                    return False
                else:
                    v, i = v ** 2, i + 1
    return True

这里使用了快速幂，否则会比较慢。

Rabin-Miller 测试是随机化算法，有一定概率会将合数误判为素数。在这里我们循环了 64 次，随机选取 64 个不同的基数进行测试，发生错误的概率小于 2^-128。

素数生成

素数筛法（生成小素数列表）

import math
def make_prime_list(n):
    assert n > 1
    prime_list = [i > 1 for i in range(n+1)]
    for i in range(2, int(math.sqrt(n))+1):
        if prime_list[i]: 
            for j in range(i * 2, n+1, i):
                prime_list[j] = False
    return prime_list

结果是长度为 n + 1 的 list。对于数 0 <= x <= n，如果 x 是素数那么 prime_list[x] 的值为 True，如果不是 prime_list[x] 的值为 False。

n 太大会导致内存不足。

生成大素数

生成随机数，然后通过 Rabin-Miller 测试判断是否是素数。

欧几里得算法（求两数最大公因子）

python 有内置的 gcd: math.gcd

def gcd(a, b):
    assert a >= 0 and b >= 0
    while a > 0:
        a, b = b % a, a
    return b

扩展欧几里得算法

一下算法出自参考书目第 112 页（原版书籍 171 页）

def extended_gcd(a, b):
    assert a >= 0 and b >= 0
    c, d = a, b
    uc, vc, ud, vd = 1, 0, 0, 1
    while c != 0:
        # 不变式 uc × a + vc × b = c ∧ ud × a + vd ∧ b = d
        # q = d // c
        # c, d = d - q * c, c
        q, c, d = divmod(d, c), c
        uc, vc, ud, vd = ud - q * uc, vd - q * vc, uc, vc
    return d, ud, vd

扩展欧几里得算法返回三个值，a b 两数的最大公约数gcd，还有 x，y。x 和 y 满足 x × a + y × b = gcd。本算法是一个递推的算法（也可以写成递归形式）。在递推循环里 c 和 d 的值不断的改变，但是这所有的 c 和 d 的最大公约数是不变的。c 和 d 的变化规律和欧几里得算法是一样的。c 每次变为 d % c，而 d 变为 c 原来的值。例如 8 和 28，最大公约数是 4。

	c	d	x	y
初始	8	28	-3	1
第一次迭代后	4	8	1	0
第二次迭代后	0	4	0	1

上表中的三个 x 分别设为 x_0，x_1，x_2。其他也一样那么

c_0 × x_0 + d_0 × y_0 = c_1 × x_1 + d_1 × y_1 = c_2 × x_2+ d_2 × y_2 = gcd = 4。

这个算法可以写的更加简洁如下：

def extended_gcd(a, b):
    assert a >= 0 and b >= 0
    ua, va, ub, vb = 1, 0, 0, 1
    while a != 0:
        (q, a), b = divmod(b, a), a
        ua, va, ub, vb = ub - q * ua, vb - q * va, ua, va
    return b, ub, vb

参考资料：

Niels Ferguson, Bruce Schneier, Tadayoshi Kohno. 《密码工程：原理与应用》,机械工业出版社; 第1版 (2018年1月1日),ISBN:9787111574354