HTB University CTF 2023 - RSA with msb of CRT exponents via Coppersmith

htb23_mayday.sage

from Crypto.Util.number import getPrime, GCD, bytes_to_long, long_to_bytes
from sage.all import PolynomialRing, RationalField, Integers, inverse_mod, Zmod, gcd

N = 0x78fb80151a498704541b888b9ca21b9f159a45069b99b04befcb0e0403178dc243a66492771f057b28262332caecc673a2c68fd63e7c850dc534a74c705f865841c0b5af1e0791b8b5cc55ad3b04e25f20dedc15c36db46c328a61f3a10872d47d9426584f410fde4c8c2ebfaccc8d6a6bd1c067e5e8d8f107b56bf86ac06cd8a20661af832019de6e00ae6be24a946fe229476541b04b9a808375739681efd1888e44d41196e396af66f91f992383955f5faef0fc1fc7b5175135ab3ed62867a84843c49bdf83d0497b255e35432b332705cd09f01670815ce167aa35f7a454f8b26b6d6fd9a0006194ad2f8f33160c13c08c81fe8f74e13e84e9cdf6566d2f
e = 0x4b3393c9fe2e50e0c76920e1f34e0c86417f9a9ef8b5a3fa41b381355
c = 0x17f2b5a46e4122ff819807a9d92b6225c483cf93c9804381098ecd6b81f4670e94d8930001b760f1d26bc7aa7dda48c9e12809d20b33fdb4c4dd9190b105b7dab42e932b99aaff54023873381e7387f1b2b18b355d4476b664d44c40413d82a10635fe6e7322543943aed2dcfbe49764b8da70edeb88d6f63ee47f025be5f2f38319611ab74cd5db6f90f60870ecbb57a884f821d873db06aadf0e61ff74cc7d4c8fc1e527dba9b205220c6707f750822c675c530f8ad6956e41ab80911da49c3d6a7d27e93c44ba5968f2f47a9c5a2694c9d6da245ceffe9cab66b6043774f446b1b08ee4739d3cc716b87c8225a84d3c4ea2fdf68143d09f062c880a870554
dp = 0x59a2219560ee56e7c35f310a4d101061aa61e0ae4eae7605eb63784209ad488b4ed161e780811edd61bf593e2d385beccfd255b459382d8a9029943781b540e7
dq = 0x39719131fbfd8afbc972ca005a430d080775bf1a5b3e8b789aba5c5110a31bd155ff13fba1019bb6cb7db887685e34ca7966a891bfad029b55b92c11201559e5

i = 512
dpM = dp
dqM = dq
A = (pow(2, 2*i) * pow(e, 2) * dpM * dqM)//N + 1


x = PolynomialRing(RationalField(), 'x').gen()
C = (1 - A*(N-1)) % e
f = x**2 - C*x + A
roots = f.roots()
if roots == []:
    f = x**2 - (C + e)*x + A
    roots = f.roots()
k = roots[0][0]
l = roots[1][0]
assert k*l == A


x = PolynomialRing(Zmod(k*N), 'x', implementation='NTL').gen()
# Assume k is the coefficient of dp
a_small = (e * dpM * (2**i) + k - 1) * inverse_mod(e, k*N)
a_small = int(a_small)
f = x + a_small
my_dpL = f.small_roots(X=2**i-1, beta=0.5)

# Is actually with dq
if my_dpL == []:
    a_small = (e * dqM * (2**i) + k - 1) * inverse_mod(e, k*N)
    a_small = int(a_small)
    f = x + a_small
    my_dqL = f.small_roots(X=2**i-1, beta=0.5)[0]
    p = gcd(f.subs(x=my_dqL), N)
    print(f'{p = }')

else:
    my_dpL = my_dpL[0]
    p = gcd(f.subs(x=my_dpL), N)
    print(f'{p = }')

p = int(p)
assert N % p == 0
q = N // p
phi = (p-1) * (q-1)
d = inverse_mod(e, phi)
print(long_to_bytes(pow(c, d, N)))