-
Notifications
You must be signed in to change notification settings - Fork 14
/
cm_factor.sage
93 lines (70 loc) · 3.43 KB
/
cm_factor.sage
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
89
90
91
92
93
import sys
from utils import *
sys.setrecursionlimit(100000000)
def cm_factor(N, D, B=32, debug=True):
"""
Implements a simplified version of Cheng's 4p - 1 elliptic curve complex multiplication based factorization algorithm.
Targets moduli where one (or many) prime factors satisfies the equality 4p_i - 1 = D * s^2, where D is a squarefree integer
and s is a randomly generated one between an upper and lower bound.
Original paper title: I want to break square-free: The 4p - 1 factorization method and its RSA backdoor viability
Link: https://crocs.fi.muni.cz/_media/public/papers/2019-secrypt-sedlacek.pdf
:param N: integer to be factored
:param D: squarefree integer satisfying 4p - 1 = D * s^2
:param B: number of iterations to run the algorithm for
:param debug: switches debugging information on/off
:return: a tuple corresponding to p, q, or failure (-1)
"""
# If D is not squarefree then we terminate immediately.
assert D.is_squarefree(), "D must be squarefree."
# Computes the -Dth Hilbert polynomial modulo N and quotient ring Q = Z_N[x] / <H_{-D, N}>
Z_N = Zmod(N)
P = Z_N[x]
H = P(hilbert_class_polynomial(-D))
Q = P.quotient_ring(H)
# j is the equivalence class corresponding to [X] in Q.
j = Q(x)
# The paper claims that we can treat both 1728 - j and H as polynomials in Z_N[X] and calculate the inverse of
# 1728 - j and H via egcd. This doesn't quite work off the shelves, so we instead accomplish this by treating
# 1728 - j as an element of Q = Z_N[x] / <H_{-D, n}> and lifting it back into the base ring.
# Sage implements this via the .lift() method.
if debug:
print("Calculating inverse mod of 1728 - j with H...")
try:
k = j * polynomial_inv_mod((1728 - j).lift(), H)
except ValueError as e:
r = gcd(int(e.args[1].lc()), N)
return int(r), int(N // r)
# Constructs an elliptic curve described by the equation y^2 = x^3 + ax + b where a = 3k and b = 2k over Q.
# This is so we can calculate the division polynomial psi_n(a, b, x_i) later.
E = EllipticCurve(Q, [3*k, 2*k])
for i in range(B):
# Obtains a random element from Z_N for the division polynomial.
x_i = Z_N.random_element()
if debug:
print("Calculating division polynomial with x_i...")
# Calculates the division polynomial modulo n: psi_n(a, b, x_i)
z = E.division_polynomial(N, x=Q(x_i))
# If our egcd implementation throws an exception, then we know that the r polynomial during the computation
# has no inverse; we return this polynomial in our exception and then take the gcd of its leading coefficient with N.
# Otherwise we take the gcd of d and N normally.
try:
d, _, _ = polynomial_egcd(z.lift(), H)
except ValueError as e:
r = gcd(int(e.args[1].lc()), N)
return int(r), int(N // r)
r = gcd(int(d), N)
# Ensure that the factorization is nontrivial: i.e. r != 1 and r != N
if 1 < r < N:
return r, N // r
return -1
if __name__ == "__main__":
D = 419
# Generates a modulus N = pq where 4p - 1 = Ds^2.
print("Generating modulus...")
p = generate_cm_prime(D)
q = random_prime(1<<1024, proof=False)
N = p * q
print("Factoring...")
r, s = cm_factor(N, D)
assert r * s == N
print(f"Factorization successful! N = {r} * {s}")