This was an interesting RSA with weak prime generation challenge. The main point of the challenge is to factor the RSA modulus, which is constructed by multiplying two primes that differ only by a random number of bytes in the middle of the number.
Here is the challenge code:
import gmpy2
import random
from Crypto.Util.number import bytes_to_long
def getPrimes(size):
half = random.randint(16, size // 2 - 8)
rand = random.randint(8, half - 1)
sizes = [rand, half, size - half - rand]
while True:
p, q = 0, 0
for s in sizes:
p <<= s
q <<= s
chunk = random.getrandbits(s)
p += chunk
if s == sizes[1]:
chunk = random.getrandbits(s)
q += chunk
p |= 2**(size - 1) | 2**(size - 2) | 1
q |= 2**(size - 1) | 2**(size - 2) | 1
if gmpy2.is_prime(p) and gmpy2.is_prime(q):
return p, q
e = 65537
p, q = getPrimes(1024)
N = p*q
phi = (p-1)*(q-1)
d = gmpy2.invert(e, phi)
m = bytes_to_long(open("flag.txt", "rb").read())
c = pow(m, e, N)
print("N =", hex(N))
print("c =", hex(c))
So essentially we can parameterize the primes using two numbers, the top $l$ bits and the bottom $k$ bits. For the purpose of our math, we represent the two primes as follows: \begin{align} p & = 2^{1024-l} a + 2^k x + b \\ q & = 2^{1024-l} a + 2^k y + b \\ n & = p q \end{align}
Solving for $a$ is easy, since $n$ is almost a square the top $l$ bits of $\sqrt{n}$ is a going to be $a$ plus or minus a small margin of error. Solving for $b$ is a bit more complicated, but it follows the same principle except instead of taking the square root of $n$ over the integers we need to solve for
\begin{align} b^2 & \equiv n \mod 2^{512} \end{align}This follows directly from the definition of $p$ and $q$ above and the fact that $l < k$
\begin{align} p & \equiv b \mod 2^{k} \\ q & \equiv b \mod 2^{k} \\ p q & \equiv b^2 \mod 2^{k} \end{align}As it turns out there is a simple and easy to implement algorithm to find the square roots of a number modulo a power of two. Essentially what we do is solve for the simple case of $\mod 2$ and then repeatedly Hensel Lift the solution until we have a set of solutions $\mod 2^{k}$. Here is the code:
import gmpy2
from Crypto.Util.number import long_to_bytes
def find_square_roots_mod_2power(n, power) :
roots = list([1])
mexp = 1
while mexp <= power :
old_m = 2**(mexp-1)
m = 2**mexp
target = n%m
new_roots = list()
for root in roots :
root1 = root
if (root1*root1) % m == target :
new_roots.append(root1)
root2 = root + old_m
if (root2*root2) % m == target :
new_roots.append(root2)
roots = new_roots
mexp += 1
return roots
Now that we have a way of solving for $a$ and $b$, we need to glue it all together and use this knowledge to factor $n$. We can do this by solving for $x+y$ instead of solving for $x$ and $y$ separately which turns it into a single-variable equation.
For convenience, now that we know $a$ and $b$ lets replace them with a new variable $c$ so that $p$ and $q$ are now defined as \begin{align} p & = 2^k x + c \\ q & = 2^k y + c \\ n &= 2^{2k} x y + 2^{k} (x + y) c + c^2 \end{align}
Using the knowledge that $x$ and $y$ are both $k$ bits, we can do the old trick of removing the term with the large exponent by taking the modulus and ending up with an equation that we know to be valid over the integers. Specifically,
\begin{align} r & = (n - c^2) / 2^{k} \\ r & = 2^k x y + (x + y) \\ \\ r & \equiv c(x + y) & \mod 2^k \\ c^{-1} r & \equiv x + y &\mod 2^k \end{align}Which gives us the Fermat factorization of $n$ and thus its two prime factors.
Also, it is worth mentioning that we don't actually know $k$ and $l$, which isn't a problem since we can easily brute force them.
import gmpy2
def factor_not_enough_lsb(n) :
for k in range(300,520) :
for l in range(1,k-1) :
a = int(gmpy2.isqrt(n) >> (1024-l))
for b in find_square_roots_mod_2power(n, k) :
c = a * 2**(1024-l) + b
r = (n-c*c) // (2**k)
x_plus_y = int(r*gmpy2.powmod(c, -1, 2**k)) % 2**k
fermatSquareA = (x_plus_y*(2**k) + 2*c) // 2
if fermatSquareA**2 - n <= 0 :
continue
p = gmpy2.isqrt(fermatSquareA**2 - n) + fermatSquareA
if n % p != 0 :
continue
q = n // p
return p, q
print('All attempts failed :(')
return None
N = 0xf7e6ddd1f49d9f875904d1280c675e0e03f4a02e2bec6ca62a2819d521441b727d313ec1b5f855e3f2a69516a3fea4e953435cbf7ac64062dd97157b6342912b7667b190fad36d54e378fede2a7a6d4cc801f1bc93159e405dd6d496bf6a11f04fdc04b842a84238cc3f677af67fa307b2b064a06d60f5a3d440c4cfffa4189e843602ba6f440a70668e071a9f18badffb11ed5cdfa6b49413cf4fa88b114f018eb0bba118f19dea2c08f4393b153bcbaf03ec86e2bab8f06e3c45acb6cd8d497062f5fdf19f73084a3a793fa20757178a546de902541dde7ff6f81de61a9e692145a793896a8726da7955dab9fc0668d3cfc55cd7a2d1d8b631f88cf5259ba1
c = 0xf177639388156bd71b533d3934016cc76342efae0739cb317eb9235cdb97ae50b1aa097f16686d0e171dccc4ec2c3747f9fbaba4794ee057964734835400194fc2ffa68a5c6250d49abb68a9e540b3d8dc515682f1cd61f46352efc8cc4a1fe1d975c66b1d53f6b5ff25fbac9fa09ef7a3d7e93e2a53f9c1bc1db30eed92a30586388cfef4d68516a8fdebe5ebf9c7b483718437fcf8693acd3118544249b6e62d30afa7def37aecf4da999c1e2b686ca9caca1b84503b8794273381b51d06d0dfb9c19125ce30e67a8cf72321ca8c50a481e4b96bbbc5b8932e8d5a32fa040c3e29ced4c8bf3541e846f832a7f9406d263a592c0a9bce88be6aed043a9867a7
p, q = factor_not_enough_lsb(N)
phi=(p-1)*(q-1)
e = 65537
N = p*q
phi = (p-1)*(q-1)
d = gmpy2.invert(e, phi)
print(long_to_bytes(pow(c, d, N)))