整数分解

发表于 分类于 阅读次数： 阅读次数：

关于整数分解算法学习记录

参考连接：https://en.wikipedia.org/wiki/Integer_factorization

1. 试除法

该方法思路是将每个数都尝试作为整数的因子去寻找，速度慢，没有什么可说的

2. Wheel factorization

基于剔除素数的倍数，筛选出素数后进行试除。

优化方法：1.埃拉托色尼筛 2.欧拉筛法

3. Pollard’s rho algorithm

参考连接：

[http://jayxv.github.io/2019/11/11/密码学学习笔记之浅析Pollard’s rho algorithm及其应用/](https://jayxv.github.io/2019/11/11/密码学学习笔记之浅析Pollard’s rho algorithm及其应用/)

原理有兴趣可以研究一下我写的另一篇博客《algorithm算法原理证明》

算法的核心思想是通过随机选择一个起始点和一个随机的多项式来寻找两个非平凡的因子。在每一步迭代中，算法计算多项式的值，并使用 Floyd 的环检测算法（Floyd’s cycle-finding algorithm）来寻找环。如果环的长度为偶数，则找到了一个非平凡的因子，否则，算法会继续尝试其他的起始点和多项式，直到找到两个因子为止。

该算法的具体实现如下：

  1. 选择一个起点x0，然后根据一个递推公式生成一系列数列：
    $$
    x_i = f(x_{i-1})
    $$
    其中f是一个特定的函数。在Pollard’s ρ算法中，函数f通常是下面这个式子：
    $$
    f(x) = x^2 + c \pmod n
    $$
    其中c是一个随机数，n是需要分解的大整数。

  2. 然后在数列中选择两个不同的位置i和j（i!=j）。

  3. 对这两个位置的元素进行求差运算，即：
    $$
    d = gcd(|xi - xj|, n)
    $$
    如果d不等于1或者n，则说明已经找到了n的一个因子，结束算法。

  4. 如果没有找到因子，则继续执行步骤1，生成新的数列，直到找到一个因子为止。

Pollard’s ρ算法是一种随机算法，因此其运行时间取决于随机数的选择。在实际应用中，该算法的平均时间复杂度约为O(sqrt(n))。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def gcd(a, b):
    while b:
        a, b = b, a%b
    return a

def mapx(x):
    x=(x*x+1)%n
    return x

def pollard_rho(x1,x2):
    while 1:
        x1 = mapx(x1)
        x2 = mapx(mapx(x2))
        p = gcd(a - i, n)
        if p != 1:
            return p
n=...

while (n!=1):
    p = pollard_rho(2,2)
    print p
    n = n / p

4.Euler’s factorization method

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
import math

def euler_factorize(n):
    if is_prime(n):
        print("Number is not factorable")
        return []
    
    for a in range(1, math.ceil(math.sqrt(n))):
        b2 = n - a*a
        b = math.floor(math.sqrt(b2))
        if b*b == b2:
            break
    else:
        print("Failed to find any expression for n as sum of squares")
        return []
    
    for c in range(a+1, math.ceil(math.sqrt(n))):
        d2 = n - c*c
        d = math.floor(math.sqrt(d2))
        if d*d == d2:
            break
    else:
        print("Failed to find a second expression for n as sum of squares")
        return []
    
    A, B = c-a, c+a
    C, D = b-d, b+d
    k = math.gcd(A, C) // 2
    h = math.gcd(B, D) // 2
    l = math.gcd(A, D) // 2
    m = math.gcd(B, C) // 2
    factor1 = k*k + h*h
    factor2 = l*l + m*m
    return [factor1, factor2]

5. Lenstra elliptic-curve factorization算法

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
import numpy as np

def is_prime(n):
    if n <= 1:
        return False
    elif n <= 3:
        return True
    elif n % 2 == 0 or n % 3 == 0:
        return False
    i = 5
    while i*i <= n:
        if n % i == 0 or n % (i+2) == 0:
            return False
        i += 6
    return True

def scalar_multiply(k, P, a, p):
    if k == 0:
        return None
    elif k == 1:
        return P
    elif k % 2 == 0:
        Q = scalar_multiply(k // 2, P, a, p)
        return add_points(Q, Q, a, p)
    else:
        Q = scalar_multiply((k-1) // 2, P, a, p)
        return add_points(add_points(Q, Q, a, p), P, a, p)

def add_points(P, Q, a, p):
    if P == Q:
        s = ((3*P[0]*P[0] + a) * pow(2*P[1], p-2, p)) % p
    else:
        s = ((Q[1] - P[1]) * pow(Q[0] - P[0], p-2, p)) % p
    x = (s*s - P[0] - Q[0]) % p
    y = (s*(P[0] - x) - P[1]) % p
    return (x, y)

def lenstra_factorize(n, t=1000):
    if n < 2:
        return []
    if n == 2:
        return [2]
    if n % 2 == 0:
        return [n//2, 2]
    if is_prime(n):
        return [n]

    for i in range(1, t+1):
        a = np.random.randint(1, n)
        x = np.random.randint(0, n)
        y = np.random.randint(0, n)
        P = (x, y)
        Q = scalar_multiply(2, P, a, n)
        d = np.gcd(Q[0] - P[0], n)
        if d != 1 and d != n:
            return [d, n//d]
    return []

print(lenstra_factorize(561))