d0_bit_length=530
d_bit_length=1024
N=0xcc5b706f373a79c680cec9527aac573fd435129cf16c23334085bf97832e5a6c78b633c2f244b12a62f87ec5295dd89fcf3c808c39e45a9afdbda2f8d2d0b50d61b685c0fe9eb41a7018a40f98892f96d738e2a4e740d4e507bcbd07f68c1ecb2ca10bd780ce65265a7e4da00f1031a5db9d038878a29a5ffefcaf2119720005
e=65537
d0=0x20142ae2802b877eb4dfa8a462e7d017c4d348181c367fd1a661ec9b6bbcca9dcb6601ccb6c10416b7f3c20129527346bbc136ee60f9945125cba03a9bba3720f7411
n = N.bit_length()
t_ = e.bit_length() - 1
alpha = t_ / n
beta = d_bit_length / n
M = 2 ** d0_bit_length
delta = (d_bit_length - d0_bit_length) / n
def _bm_6(N, e, d_bit_length, d0, d0_bit_length, M, m, t):
y, z = ZZ["y", "z"].gens()
f = y * (N - z) - e * d0 + 1
Y = e # Equivalent to N^alpha
Z = int(3 * RR(N) ** (1 / 2))
for y0, z0 in modular_bivariate(f, e * M, m, t, Y, Z):
phi = N - z0
d = pow(e, -1, phi)
print(d)
if pow(pow(2, e, N), d, N) == 2:
factors = factorize(N, phi)
if factors:
return *factors, d
return None
def factorize(N, phi):
"""
Recovers the prime factors from a modulus if Euler's totient is known.
This method only works for a modulus consisting of 2 primes!
:param N: the modulus
:param phi: Euler's totient, the order of the multiplicative group modulo N
:return: a tuple containing the prime factors, or None if the factors were not found
"""
s = N + 1 - phi
d = s ** 2 - 4 * N
p = int(s - isqrt(d)) // 2
q = int(s + isqrt(d)) // 2
return p, q if p * q == N else None
def modular_bivariate(f, eM, m, t, Y, Z, roots_method="groebner"):
"""
Computes small modular roots of a bivariate polynomial.
More information: Blomer J., May A., "New Partial Key Exposure Attacks on RSA" (Section 6)
:param f: the polynomial
:param eM: the modulus
:param m: the parameter m
:param t: the parameter t
:param Y: an approximate bound on the y roots
:param Z: an approximate bound on the z roots
:param roots_method: the method to use to find roots (default: "groebner")
:return: a generator generating small roots (tuples of y and z roots) of the polynomial
"""
f = f.change_ring(ZZ)
pr = f.parent()
y, z = pr.gens()
shifts = []
for i in range(m + 1):
for j in range(i + 1):
g = y ** j * eM ** i * f ** (m - i)
shifts.append(g)
for j in range(1, t + 1):
h = z ** j * eM ** i * f ** (m - i)
shifts.append(h)
L, monomials = create_lattice(pr, shifts, [Y, Z])
L = reduce_lattice(L)
polynomials = reconstruct_polynomials(L, f, eM ** m, monomials, [Y, Z])
for roots in find_roots(pr, polynomials, method=roots_method):
yield roots[y], roots[z]
def create_lattice(pr, shifts, bounds, order="invlex", sort_shifts_reverse=False, sort_monomials_reverse=False):
"""
Creates a lattice from a list of shift polynomials.
:param pr: the polynomial ring
:param shifts: the shifts
:param bounds: the bounds
:param order: the order to sort the shifts/monomials by
:param sort_shifts_reverse: set to true to sort the shifts in reverse order
:param sort_monomials_reverse: set to true to sort the monomials in reverse order
:return: a tuple of lattice and list of monomials
"""
if pr.ngens() > 1:
pr_ = pr.change_ring(ZZ, order=order)
shifts = [pr_(shift) for shift in shifts]
monomials = set()
for shift in shifts:
monomials.update(shift.monomials())
shifts.sort(reverse=sort_shifts_reverse)
monomials = sorted(monomials, reverse=sort_monomials_reverse)
L = matrix(ZZ, len(shifts), len(monomials))
for row, shift in enumerate(shifts):
for col, monomial in enumerate(monomials):
L[row, col] = shift.monomial_coefficient(monomial) * monomial(*bounds)
monomials = [pr(monomial) for monomial in monomials]
return L, monomials
def reduce_lattice(L, delta=0.8):
"""
Reduces a lattice basis using a lattice reduction algorithm (currently LLL).
:param L: the lattice basis
:param delta: the delta parameter for LLL (default: 0.8)
:return: the reduced basis
"""
return L.LLL(delta)
def reconstruct_polynomials(B, f, modulus, monomials, bounds, preprocess_polynomial=lambda x: x, divide_gcd=True):
"""
Reconstructs polynomials from the lattice basis in the monomials.
:param B: the lattice basis
:param f: the original polynomial (if set to None, polynomials will not be divided by f if possible)
:param modulus: the original modulus
:param monomials: the monomials
:param bounds: the bounds
:param preprocess_polynomial: a function which preprocesses a polynomial before it is added to the list (default: identity function)
:param divide_gcd: if set to True, polynomials will be pairwise divided by their gcd if possible (default: True)
:return: a list of polynomials
"""
divide_original = f is not None
modulus_bound = modulus is not None
polynomials = []
for row in range(B.nrows()):
norm_squared = 0
w = 0
polynomial = 0
for col, monomial in enumerate(monomials):
if B[row, col] == 0:
continue
norm_squared += B[row, col] ** 2
w += 1
assert B[row, col] % monomial(*bounds) == 0
polynomial += B[row, col] * monomial // monomial(*bounds)
print(661)
# Equivalent to norm >= modulus / sqrt(w)
#if modulus_bound and norm_squared * w >= modulus ** 2:
#continue
print(662)
polynomial = preprocess_polynomial(polynomial)
if divide_original and polynomial % f == 0:
polynomial //= f
if divide_gcd:
for i in range(len(polynomials)):
g = gcd(polynomial, polynomials[i])
# TODO: why are we only allowed to divide out g if it is constant?
if g != 1 and g.is_constant():
polynomial //= g
polynomials[i] //= g
if polynomial.is_constant():
continue
polynomials.append(polynomial)
return polynomials
def find_roots_univariate(x, polynomial):
"""
Returns a generator generating all roots of a univariate polynomial in an unknown.
:param x: the unknown
:param polynomial: the polynomial
:return: a generator generating dicts of (x: root) entries
"""
if polynomial.is_constant():
return
for root in polynomial.roots(multiplicities=False):
if root != 0:
yield {x: int(root)}
def find_roots_gcd(pr, polynomials):
"""
Returns a generator generating all roots of a polynomial in some unknowns.
Uses pairwise gcds to find trivial roots.
:param pr: the polynomial ring
:param polynomials: the reconstructed polynomials
:return: a generator generating dicts of (x0: x0root, x1: x1root, ...) entries
"""
if pr.ngens() != 2:
return
x, y = pr.gens()
for i in range(len(polynomials)):
for j in range(i):
g = gcd(polynomials[i], polynomials[j])
if g.degree() == 1 and g.nvariables() == 2 and g.constant_coefficient() == 0:
# g = ax + by
a = int(g.monomial_coefficient(x))
b = int(g.monomial_coefficient(y))
yield {x: b, y: a}
yield {x: -b, y: a}
def find_roots_groebner(pr, polynomials):
"""
Returns a generator generating all roots of a polynomial in some unknowns.
Uses Groebner bases to find the roots.
:param pr: the polynomial ring
:param polynomials: the reconstructed polynomials
:return: a generator generating dicts of (x0: x0root, x1: x1root, ...) entries
"""
# We need to change the ring to QQ because groebner_basis is much faster over a field.
# We also need to change the term order to lexicographic to allow for elimination.
gens = pr.gens()
s = Sequence(polynomials, pr.change_ring(QQ, order="lex"))
while len(s) > 0:
G = s.groebner_basis()
if len(G) == len(gens):
roots = {}
for polynomial in G:
vars = polynomial.variables()
if len(vars) == 1:
for root in find_roots_univariate(vars[0], polynomial.univariate_polynomial()):
roots |= root
if len(roots) == pr.ngens():
yield roots
return
G = Sequence(s, pr.change_ring(ZZ, order="lex")).groebner_basis()
vars = tuple(map(lambda x: var(x), gens))
for solution_dict in solve([polynomial(*vars) for polynomial in G], vars, solution_dict=True):
found = False
roots = {}
for i, v in enumerate(vars):
s = solution_dict[v]
if s.is_constant():
if not s.is_zero():
found = True
roots[gens[i]] = int(s) if s.is_integer() else int(s) + 1
else:
roots[gens[i]] = 0
if found:
yield roots
return
return
else:
# Remove last element (the biggest vector) and try again.
s.pop()
def find_roots_resultants(gens, polynomials):
"""
Returns a generator generating all roots of a polynomial in some unknowns.
Recursively computes resultants to find the roots.
:param polynomials: the reconstructed polynomials
:param gens: the unknowns
:return: a generator generating dicts of (x0: x0root, x1: x1root, ...) entries
"""
if len(polynomials) == 0:
return
if len(gens) == 1:
if polynomials[0].is_univariate():
yield from find_roots_univariate(gens[0], polynomials[0].univariate_polynomial())
else:
resultants = [polynomials[0].resultant(polynomials[i], gens[0]) for i in range(1, len(gens))]
for roots in find_roots_resultants(gens[1:], resultants):
for polynomial in polynomials:
polynomial = polynomial.subs(roots)
if polynomial.is_univariate():
for root in find_roots_univariate(gens[0], polynomial.univariate_polynomial()):
yield roots | root
def find_roots_variety(pr, polynomials):
"""
Returns a generator generating all roots of a polynomial in some unknowns.
Uses the Sage variety (triangular decomposition) method to find the roots.
:param pr: the polynomial ring
:param polynomials: the reconstructed polynomials
:return: a generator generating dicts of (x0: x0root, x1: x1root, ...) entries
"""
# We need to change the ring to QQ because variety requires a field.
s = Sequence([], pr.change_ring(QQ))
for polynomial in polynomials:
s.append(polynomial)
I = s.ideal()
dim = I.dimension()
if dim == -1:
s.pop()
elif dim == 0:
for roots in I.variety(ring=ZZ):
yield {k: int(v) for k, v in roots.items()}
return
def find_roots(pr, polynomials, method="groebner"):
"""
Returns a generator generating all roots of a polynomial in some unknowns.
The method used depends on the method parameter.
:param pr: the polynomial ring
:param polynomials: the reconstructed polynomials
:param method: the method to use, can be "groebner", "resultants", or "variety" (default: "groebner")
:return: a generator generating dicts of (x0: x0root, x1: x1root, ...) entries
"""
if pr.ngens() == 1:
for polynomial in polynomials:
yield from find_roots_univariate(pr.gen(), polynomial)
else:
# Always try this method because it can find roots the others can't.
yield from find_roots_gcd(pr, polynomials)
print(method)
if method == "groebner":
yield from find_roots_groebner(pr, polynomials)
elif method == "resultants":
yield from find_roots_resultants(pr.gens(), polynomials)
elif method == "variety":
yield from find_roots_variety(pr, polynomials)
if e < RR(N) ** (7 / 8) and RR(N) ** (1 / 6 + 1 / 3 * sqrt(1 + 6 * alpha)) <= M:
t = None
m=10
t = int((2 / 3 * (1 - delta - alpha) / (2 * alpha - 1))) if t is None else t
_bm_6(N, e, d_bit_length, d0, d0_bit_length, M, m, t)