Challenge:

from Crypto.Util.number import *
from secret import flag
z = 567
p = getPrime(1024)
q = getPrime(1024)
n = p*q
c = pow(bytes_to_long(flag), 65537, n)
tot = (p-1) * (q-1)
d = int(pow(65537, -1, tot))
dinv = int(pow(d, -1, n))

h = int(dinv >> z)
hpq = (int((p+q)>> z))

with open('out.txt', 'w+') as f:
    f.write(f'{n=}\n')
    f.write(f'{h=}\n')
    f.write(f'{hpq=}\n')
    f.write(f'{c=}\n')

n=13986357905153484822874300783445968480194277882812317554826224241536479785567487956712558237728345348661360577246137576216953724039680969623887884690471844396542763308129517234365819619617071449273126659007918716307793788623728052337632935762139796688014791419718949572448772521789488223910450877828732015095423443037519388747356327730350934152781671783952028215703864406564741666179193772037496984854699143314813242721157017296866888522135989818414587193505121794302821401677072507471357592358012342178011963104524959087968374300060349343826214249928530346877968114749229074874962737714935221065368318487049394644831
h=10474216468878927114435400909130676124750910912012236182806861194655854223324539867768381265996955193355030239325750528328250897464859373863289680002879536341349759323910048168674147097644573874679268018966497862685092382336865554114348153248267599439087357199554652601126191061921516650448119261614064051599968120061991607030873881013657693987836636730528537557619595799676312850875727477092697270452300532360780188724484703363561848754770976459
hpq=492124417091708682668644108145880307537308922842816506360717440112116492381514432506339907757228214359689270777951081610062506962769167209
c=4715651972688371479449666526727240348158670108161494767004202259402013317642418593561463200947908841531208327599049414587586292570298317049448560403558027904798159589477994992384199008976859139072407664659830448866472863679123027179516506312536814186903687358198847465706108667279355674105689763404474207340186200156662095468249081142604074167178023479657021133754055107459927667597604156397468414872149353231061997958301747136265344906296373544580143870450924707559398134384774201700278038470171319329716930036843839101955981274793386611943442507144153946307781795085665793554799349509983282980591388585613674226899



Solve:

The 2 hints given to us are the msb of dinv and the msb of p+q.

First I bit-shift them back.

Then denote the 2 unknowns lsbs x and y, so:

\[\text{dinv} = \text{h} + \text{x}\] \[\text{p+q} = \text{hpq} + \text{y}\]

Initial equation:

\[\text{d} \cdot \text{dinv} - 1 \equiv 0 \text{ (mod n)}\]

An equation for d (some unknown k where 0<k<e):

\[\text{d} = \frac{1 + \text{k} \cdot \text{phi}}{\text{e}}\]

Equation for phi:

\[\text{phi} = \text{n} + 1 - (\text{p}+\text{q}) = \text{n} + 1 - (\text{hpq} + \text{y})\]

Sub phi into d equation:

\[\text{d} = \frac{1 + \text{k} \cdot (\text{n} + 1 - (\text{hpq} + \text{y}))}{\text{e}}\]

Now sub d and dinv into initial equation:

\[\big(\frac{1 + \text{k} \cdot (\text{n} + 1 - (\text{hpq} + \text{y}))}{\text{e}} \big) \cdot (\text{h} + \text{x}) - 1 \equiv 0 \text{ (mod n)}\] \[\big(1 + \text{k} \cdot (\text{n} + 1 - (\text{hpq} + \text{y})) \big) \cdot (\text{h} + \text{x}) - \text{e} \equiv 0 \text{ (mod n)}\]

Now x and y are solvable with bivariate coppersmith, although you have to bruteforce k…

import itertools
from tqdm import trange

def bivariate(f, bounds, m=1, d=None):
    if d is None:
        d = f.degree()
    R = f.base_ring()
    N = R.cardinality()
    f_ = (f // f.lc()).change_ring(ZZ)
    f = f.change_ring(ZZ)
    l = f.lm()
    M = []
    for k in range(m+1):
        M_k = set()
        T = set((f^(m-k)).monomials())
        for mon in (f^m).monomials():
            if mon//l^k in T: 
                for extra in itertools.product(range(d), repeat=f.nvariables()):
                    g = mon * prod(map(power, f.variables(), extra))
                    M_k.add(g)
        M.append(M_k)
    M.append(set())
    shifts = Sequence([], f.parent())
    for k in range(m+1):
        for mon in M[k] - M[k+1]:
            g = mon//l^k * f_^k * N^(m-k)
            shifts.append(g)
    B, monomials = shifts.coefficients_monomials()
    factors = [monomial(*bounds) for monomial in monomials]
    for i, factor in enumerate(factors):
        B.rescale_col(i, factor)
    B = B.dense_matrix().LLL()
    B = B.change_ring(QQ)
    for i, factor in enumerate(factors):
        B.rescale_col(i, 1/factor)
    B = B.change_ring(ZZ)
    H = [h for h in B * monomials if not h.is_zero()]
    H = H[:2] # speedup
    for f1, f2 in itertools.combinations(H, r=2):
        x, y = f.parent().gens()
        x = f1.parent()(x)
        y = f1.parent()(y)
        res = f1.resultant(f2,y).univariate_polynomial()
        if res == 0:
            continue
        rs = res.roots()
        if rs:
            x = rs[0][0]
            if x != 0:
                y = f1.subs(x=x).univariate_polynomial().roots()
                if y:
                    y = y[0][0]
                    return (x, y)

n=13986357905153484822874300783445968480194277882812317554826224241536479785567487956712558237728345348661360577246137576216953724039680969623887884690471844396542763308129517234365819619617071449273126659007918716307793788623728052337632935762139796688014791419718949572448772521789488223910450877828732015095423443037519388747356327730350934152781671783952028215703864406564741666179193772037496984854699143314813242721157017296866888522135989818414587193505121794302821401677072507471357592358012342178011963104524959087968374300060349343826214249928530346877968114749229074874962737714935221065368318487049394644831
h=10474216468878927114435400909130676124750910912012236182806861194655854223324539867768381265996955193355030239325750528328250897464859373863289680002879536341349759323910048168674147097644573874679268018966497862685092382336865554114348153248267599439087357199554652601126191061921516650448119261614064051599968120061991607030873881013657693987836636730528537557619595799676312850875727477092697270452300532360780188724484703363561848754770976459
hpq=492124417091708682668644108145880307537308922842816506360717440112116492381514432506339907757228214359689270777951081610062506962769167209
c=4715651972688371479449666526727240348158670108161494767004202259402013317642418593561463200947908841531208327599049414587586292570298317049448560403558027904798159589477994992384199008976859139072407664659830448866472863679123027179516506312536814186903687358198847465706108667279355674105689763404474207340186200156662095468249081142604074167178023479657021133754055107459927667597604156397468414872149353231061997958301747136265344906296373544580143870450924707559398134384774201700278038470171319329716930036843839101955981274793386611943442507144153946307781795085665793554799349509983282980591388585613674226899
e = 65537
z = 567
h <<= z
hpq <<= z

PR.<x, y> = PolynomialRing(Zmod(n))
#for k in trange(e): 
for k in trange(8500, e): # skip
    f = (1 + k * (n + 1 - (hpq + y))) * (h + x) - e 
    roots = bivariate(f, bounds=(2**z, 2**z), m=2, d=1)
    if roots is None:
        continue
    d = pow(h + roots[0], -1, n)
    flag = int(pow(c, d, n))
    break

print(f'{k = }')
print(bytes.fromhex(f'{flag:x}'))



ETA: trivariate is possible:

load('https://raw.githubusercontent.com/Connor-McCartney/coppersmith/refs/heads/main/coppersmith.sage')

n=13986357905153484822874300783445968480194277882812317554826224241536479785567487956712558237728345348661360577246137576216953724039680969623887884690471844396542763308129517234365819619617071449273126659007918716307793788623728052337632935762139796688014791419718949572448772521789488223910450877828732015095423443037519388747356327730350934152781671783952028215703864406564741666179193772037496984854699143314813242721157017296866888522135989818414587193505121794302821401677072507471357592358012342178011963104524959087968374300060349343826214249928530346877968114749229074874962737714935221065368318487049394644831
h=10474216468878927114435400909130676124750910912012236182806861194655854223324539867768381265996955193355030239325750528328250897464859373863289680002879536341349759323910048168674147097644573874679268018966497862685092382336865554114348153248267599439087357199554652601126191061921516650448119261614064051599968120061991607030873881013657693987836636730528537557619595799676312850875727477092697270452300532360780188724484703363561848754770976459
hpq=492124417091708682668644108145880307537308922842816506360717440112116492381514432506339907757228214359689270777951081610062506962769167209
c=4715651972688371479449666526727240348158670108161494767004202259402013317642418593561463200947908841531208327599049414587586292570298317049448560403558027904798159589477994992384199008976859139072407664659830448866472863679123027179516506312536814186903687358198847465706108667279355674105689763404474207340186200156662095468249081142604074167178023479657021133754055107459927667597604156397468414872149353231061997958301747136265344906296373544580143870450924707559398134384774201700278038470171319329716930036843839101955981274793386611943442507144153946307781795085665793554799349509983282980591388585613674226899
e = 65537
z = 567
h <<= z
hpq <<= z

PR.<k, x, y> = PolynomialRing(Zmod(n))
f = (1 + k * (n + 1 - (hpq + y))) * (h + x) - e 
k, x, y = multivariate(f, (e, 2**z, 2**z), "shift_polynomials", "jacobian", m=2, d=1)[0]
d = pow(h+x, -1, n)
flag = int(pow(c, d, n))
print(bytes.fromhex(f'{flag:x}'))