肝了一天总算把大数质因数分解搞定了，这篇文章主要涉及了 Pollard rho 算法和试除法
我用了若干个质数的平方来对比这两个算法的性能，发现：
- 7e5 以上的数用 Pollard rho 算法更快，分解多大的数都不是问题
- 7e5 以下的数用试除法更快
所以最终的质因数分解是由这两个算法构成的，以下是基本的配备
from random import randint
from math import gcd, isqrt
def quick_power(basic, time, mod):
''' 二分快速幂
mod: 求余数'''
result = 1
while time:
if time & 1:
result = result * basic % mod
basic = basic ** 2 % mod
time >>= 1
return result
主函数的思路是：
- 当 n > 7e5 时，使用 Miller Rabin 算法判断 n 是不是质数，如果不是：
- 使用 Pollard rho 算法求解任意因数 a，并求出 a 对 n 的整除次数（n 更新为整除后的值）
- 递归分解出 a 的质因数；如果 n > 1，则递归分解出当前 n 的质因数
- 当 n <= 7e5 时，使用试除法进行质因数分解
至此涉及到的算法有：试除法，Miller Rabin 素性测试，Pollard rho 因数分解
Miller Rabin 素性测试对于素数 p 和随机数 a，如果两者互素，根据费马小定理有：
a 通常取预设的值，如：[2, 325, 9375, 28178, 450775, 9780504, 1795265022]
这时，如果 p 满足上式则很有可能是素数，否则一定是合数
因为素数都是奇数，所以 p - 1 一定是偶数，所以有：
可被表示成若干个因子的累乘，如果其中有因子可以被 p 整除，则 成立
即当满足以下任意式子时，p 通过 a 的素性测试：
......
为了提高准确率，通常取不同的 a，进行 5 ~ 10 次测试
def miller_rabin(p):
''' 素性测试'''
# 特判 4
if p <= 4: return p in (2, 3)
# 对 p-1 进行分解
pow_2, tmp = 0, p - 1
while tmp % 2 == 0:
tmp //= 2
pow_2 += 1
# 进行多次素性测试
for a in [2, 325, 9375, 28178, 450775, 9780504, 1795265022]:
basic = quick_power(a, tmp, p)
# a^m 是 p 的倍数或者满足条件
if basic in [0, 1, p - 1]:
# 通过 a 的素性测试
continue
# 进行 r-1 次平方
pass_flag = False
for _ in range(1, pow_2):
basic = basic ** 2 % p
# 怎样平方都是 1
if basic == 1: return False
# 通过 a 的素性测试
if basic == p - 1:
pass_flag = True
break
# 未通过 a 的素性测试
if not pass_flag: return False
# 通过所有 a 的素性测试
return True
Pollard rho 因数分解
我比较菜，就不讲解原理了，具体见这篇文章：Link ~
def pollard_rho(n):
''' 返回因数'''
if n == 4: return 2
bias = randint(3, n - 1)
# 更新函数
update = lambda i: (i ** 2 + bias) % n
x = randint(0, n - 1)
y = update(x)
# 查找序列环
while x != y:
factor = gcd(abs(x - y), n)
if factor != 1: return factor
x = update(x)
y = update(update(y))
return n
质因数分解
使用继承 dict 字典的类，主要是为了记录因数及其幂次，这个类主要有几个函数：
- add：更新因数的幂次
- count：试除并求出幂次
- main：主函数，内含质因数分解的主流程
- try_devide：试除法分解质因数
这个类可以当作函数来用，prime_factor(n) 的返回值就是一个字典
class prime_factor(dict):
''' 质因数分解
require: quick_power, miller_rabin, pollard_rho'''
def __init__(self, n):
super(prime_factor, self).__init__()
self.main(n)
def add(self, n, cnt):
# 更新因数表
self[n] = self.get(n, 0) + cnt
def count(self, n, fac):
# 试除并记录幂次
cnt = 1
n //= fac
while n % fac == 0:
cnt += 1
n //= fac
return n, cnt
def main(self, n, gain=1):
if n > 7e5:
# 米勒罗宾判素
if miller_rabin(n):
self.add(n, gain)
else:
# pollard rho 求解因数
fac = pollard_rho(n)
# 求解幂次
n, cnt = self.count(n, fac)
# 递归求解因数的因数
self.main(fac, gain=cnt)
# 递归求解剩余部分
if n > 1: self.main(n)
# 试除法求解
else:
self.try_divide(n, gain=gain)
def try_divide(self, n, gain=1):
''' 试除法分解'''
i, bound = 2, isqrt(n)
while i <= bound:
if n % i == 0:
# 计数 + 整除
n, cnt = self.count(n, i)
# 记录幂次, 更新边界
self.add(i, cnt * gain)
bound = isqrt(n)
i += 1
if n > 1: self.add(n, gain)
