0%

Python大数因数分解

每个合数都可以写成几个质数相乘的形式,其中每个质数都是这个合数的因数,叫做这个合数的分解质因数。本篇文章会用三种方法来进行因数分解。
如:$7804725584345565904628551916716033 = 19 \cdot 73 \cdot 307 \cdot 465841 \cdot 31865908033 \cdot 1234749313729$。
而$19$、$73$、$307$、$465841$、$31865908033$、 $1234749313729$均为质数。


短除法-1

从i为2开始枚举,一直枚举到$\sqrt{n}$,一旦n % i == 0成立,则i为n的因子,然后进行n //= i使运行速度加快并使i为质数才可能有n % i == 0

  • 复杂度:$O(\sqrt{n})$
  • 适用范围:$n\in[2, 10^{16}]$
1
2
3
4
5
6
7
8
9
10
11
12
13
14
def factorization(n):
i = 2
ret = []
while i * i <= n:
while n % i == 0:
ret.append(i)
n //= i
i += 1
if n > 1:
ret.append(n)
return ret

if __name__ == '__main__':
print(factorization(int(input())))

由于这种算法会试除很多合数,而我们知道合数是不可能满足n % i == 0,而在自然数中合数比质数多很多。在进行多个数的质因数分解时,这个算法就会显得更加鸡肋。

短除法-2

我们可以打表出2到$\sqrt{n}$之间的质数再进行试除,这样在解决多个数的质因数分解时才会免除大部分合数的影响。用素数筛进行打表复杂度为O(n),我们也只需要从2打表到$\sqrt{n}$即可。

  • 复杂度:$O(\sqrt{n})$
  • 适用范围:$n\in[2, 10^{12}]$
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
pri = []
MX = int(1e6)
isprime = [True] * MX
def init():
global a, MX
for i in range(2, MX):
if isprime[i]:
pri.append(i)
for j in range(i + i, MX, i):
isprime[j] = False


def factorization(n):
global pri
ret = []
for i in pri:
if i * i > n:
break
while n % i == 0:
ret.append(i)
n //= i
return ret


if __name__ == '__main__':
init()
print(factorization(int(input())))

用Miller-Rabin素性测试和离散对数Pollard_rho算法进行大数因数分解:

1975年,John M. Pollard提出。该算法时间复杂度为O($n^{\frac{1}{4}}$)。

参考资料:
Miller–Rabin primality test
Pollard’s rho algorithm

  • 复杂度:$O(n^{\frac{1}{4}})$
  • 适用范围:$n\in[2, 10^{33}]$
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
import random
from math import log, log10
from collections import Counter

def gcd(x, y):
return x if y == 0 else gcd(y, x % y)

def fpow(a, x, n):
ans = 1
while x > 0:
if x & 1:
ans = ans * a % n
a = a * a % n
x >>= 1
return ans

# there change the times of Rabin-Miller
TIMES = 10
def is_prime(n):
def check(a, n, x, t):
ret = fpow(a, x, n)
last = ret
for i in range(0, t):
ret = ret * ret % n
if ret == 1 and last != 1 and last != n - 1:
return True
last = ret
if ret != 1:
return True
return False

if not isinstance(n, int):
raise TypeError(str(n) + ' is not an integer!')
if n <= 0:
raise ValueError('%d <= 0' % n)
if n in {2, 3, 5, 7, 11}:
return True
for i in {2, 3, 5, 7, 11}:
if n % i == 0:
return False
x = n - 1
t = 0
while not x & 1:
x >>= 1
t += 1
for i in range(0, TIMES):
a = random.randint(1, n - 2)
if check(a, n, x, t):
return False
return True

def pollard_rho_2(n, c):
x = random.randint(0, n)
i, k, y = 1, 2, x
while True:
i += 1
x = (x * x) % n + c
d = gcd(y - x, n)
if d != 1 and d != n:
return d
if y == x:
return n
if i == k:
y = x
k <<= 1

def pollard_rho_1(n):
if not isinstance(n, int):
raise TypeError(str(n) + ' is not an integer!')
if n == 1:
return None
if is_prime(n):
return [n]
ans = []
p = n
while p >= n:
p = pollard_rho_2(p, random.randint(1, n - 1))
ans.extend(pollard_rho_1(p))
ans.extend(pollard_rho_1(n // p))
return ans

def factorization(n):
return Counter(pollard_rho_1(n))

if __name__ == '__main__':
n = int(input())
print('len:', len(str(n)))
print(factorization(n))

缺点是代码很长。
而且如果要用C++进行改写则记得考虑溢出(快速乘或套用大数类)。

回到开头