This challenge is called YAFM aka. "Yet Another Factorization Method.". The source code outlines a server that generates an RSA public key upon demand and then provides the encrypted flag upon request. Standard RSA with PKCS1_OAEP padding and $e = 65537$, so it seems on the surface like the right path to follow is to attack the special generate_prime function which is used to generate the two prime factors of the modulus.
def generate_prime(prime_len):
bits_len = 180
while True:
bits = random.getrandbits(bits_len)
idxs = random.sample(list(range(1, prime_len-2)), bits_len)
p = 1 | 2**(prime_len - 1) | 2**(prime_len - 2)
for i in range(bits_len):
p += (bits >> i & 1)*2**idxs[i]
if isPrime(p):
return p
This funky prime generation logic results in primes which are almost entirely full of zero bits. It starts with the number $2^{1023} + 2^{1022} + 1$, selects 180 random bits and randomly flips some of those bits to zero. In the end the expected outcome is a prime that has only 93 non-zero bits.
My approach was to use Hensel lifting to solve the polynomial $p*q = n \mod 2^{1024}$. Normally this approach would fail miserably because the polynomial actually has $2^{1024}$ solutions - but we have some additional information that we are able to exploit to trim down the number of solutions.
\begin{equation} S = \{(p, q) : p,q \in \mathbb{Z}_{1024}, pq \equiv n \mod 2^{1024}\} \end{equation}Specifically we only care about the elements $(p, q) \in S$ for which both $p$ and $q$ have a small number of bits. Assuming they are both essentially random elements of $\mathbb{Z}_{1024}$, the likelihood of one of them having 93 or less bits small number of bits is $r=\frac{\sum_{i=0}^{93}{1024 \choose i}}{2^{1024}}$. This probability is not good enough itself, the real weakness comes from the fact that the modulus we are attacking has this property for both of it primes. Thus the probability of both factors having 93 or less bits is the much smaller $r^{2}$.
So you have a clear picture of the actual numbers: \begin{align} 2^{1024}r &\approx 1.17 \times 10^{134}\\ 2^{1024}r^{2} &\approx 7.63 \times 10^{-41} \end{align}
The tricky part is that the Hensel lift will start at a small modulus and have to work its way up one step at a time to the final modulus of $2^{1024}$. Even though one of the primes may have only 93 bits it is possible that they are clustered near the low bits which could cause the prime to be trimmed out while the modulus is still small.
The way I got around this was to use a heuristic where at each step in the Hensel lift we keep at most the 1 million solutions with the smallest number of bits. This struck a good balance between performance and probability of success, and this half optimized single-threaded implementation was able to factor the generated public keys with high enough likelihood. Instead of spending time optimizing further or adding parallelism to the algorithm, I achieved parallelism by having multiple instances of the algorithm running on different public keys at the same time which did the trick.
CTF{l0w_entr0py_1s_alw4ys_4_n1ghtmar3_I_h4v3_sp0ken}
import sys
import datetime
from array import *
start = datetime.datetime.now()
n = 0x900048840fc8869113799f83a4c2c83615504969a4ac2f4fad2823cd061a42a6feed695ecbc7158bc16e8e3c3d278a99e54cbcd4abb347e942a269e16e2665347bebf3050e69a2f33380af333a9fc9c234f30c0a9e72e813a08b5ec0109ff7456edd67989281c521ff734567cec594b722bb00ab2a5ce788b4ac3a16f504dbcb9813ae34e41ffdd82e92fd331ab816bff0a5f5f57922a3267ce721b99fe7e8b558d63badcbccb486ea36f4421029ae3e6291d63651ccb3183fa422a82a7291998f595f5ccc0546c7aeba5bfe40825c2aeb2693653e1bb2affbd73b1676396830838d2ebcfde8a16ddff41eea73871fecc862798d79bbcbd8c55ed76540809811
results = list([(1,1, 1, 1)])
k = 2
def add_result(a, b, res) :
res.add( (a,b) )
def process(a, aOnes, b, bOnes, m, target, ones_counter, new_results) :
if (a*b) & m == target :
ones_counter[max(aOnes,bOnes)] += 1
new_results.append( (a, aOnes, b, bOnes) )
def advance_result(results, k, iters, carryover) :
for _ in range(iters) :
k += 1
last = 2**(k-1)
m = 2**k - 1
target = n & m
new_results = list()
ones_counter = array('i', [0 for _ in range(1025)])
for x, xOnes, y, yOnes in results :
largerX = x|last
largerY = y|last
process(x, xOnes, y, yOnes, m, target, ones_counter, new_results)
process(largerX, xOnes+1, y, yOnes, m, target, ones_counter, new_results)
process(x, xOnes, largerY, yOnes+1, m, target, ones_counter, new_results)
process(largerX, xOnes+1, largerY, yOnes+1, m, target, ones_counter, new_results)
s = 0
upper_bound = 0
while upper_bound < 1025 :
s += ones_counter[upper_bound]
if s > carryover :
break
else :
upper_bound += 1
results = list()
for entry in new_results :
x, xOnes, y, yOnes = entry
if xOnes < upper_bound and yOnes < upper_bound :
results.append(entry)
return results, k
working_set_target = 1000000
for i in range(1023) :
results, k = advance_result(results, k, 1, working_set_target)
found = False
for p, pOnes, q, qOnes in results :
if p*q == n :
print('Factors found!')
print('p =',p)
print('q =',q)
found = True
break
if not found :
print('No factors found :(')