defund/coppersmith + msolve

coppersmith.sage

import itertools

# curl -LJO 'https://msolve.lip6.fr/downloads/msolve-v0.2.4.tar.gz'
load("msolve-v0.2.4/interfaces/msolve-to-sage-file-interface.sage")

def msolve_wrap(eqs):
    return MSolveRealRoots(eqs, mspath="msolve-v0.2.4/binary/msolve", p=0)

def all_zz(xs):
    try:
        return [ZZ(x) for x in xs]
    except:
        pass


def solve_integer_equations(eqs):
    roots = msolve_wrap(eqs)
    rs = []
    for root in roots:
        r = all_zz(root)
        if r:
            rs.append(r)
    return rs


def small_roots(f, bounds, m=1, d=None):
    if not d:
        d = f.degree()

    R = f.base_ring()
    N = R.cardinality()

    f /= f.coefficients().pop(0)
    f = f.change_ring(ZZ)

    G = Sequence([], f.parent())
    for i in range(m + 1):
        base = N ^ (m - i) * f ^ i
        for shifts in itertools.product(range(d), repeat=f.nvariables()):
            g = base * prod(map(power, f.variables(), shifts))
            G.append(g)

    B, monomials = G.coefficient_matrix()
    monomials = vector(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)

    H = Sequence([], f.parent().change_ring(QQ))
    for h in filter(None, B * monomials):
        H.append(h)
        I = H.ideal()
        if I.dimension() == -1:
            H.pop()
        elif I.dimension() == 0:
            return solve_integer_equations(list(H))

    return []


if __name__ == "__main__":
    # pbctf Special Gift Revenge
    N = 123463519828344660835965296108959625188149729700517379543746606603601816029557213728343115758280318474617032830851553509268562367217512005079977122560679743955588214135519642513042848616372204042776892196887455692479457740367547908255044784496969010537283159300508751036032559594474145098337531029291955103059
    e = 85803665824396212221464259773478155183477895540333642019501498374139506738444521180470104195883386495607712971252463223185914391456070458788554837326327618859712794129800329295751565279950274474800740076285111503780662397876663144946831503522281710586712396810593754749589799811545251575782431569881989690861
    gift = 46710143823773072238724337855139753113453277386728402328859555407710009799097841900723288768522450009531777773692804519189753306306645410280934372812
    d0 = gift << 120
    R = Integers(e)
    P.<x, s> = PolynomialRing(ZZ)
    bounds = (2 ^ 120, 2 ^ 512)
    f = 1 + 2 * (x + e * d0 // N) * ((N + 1) // 2 - s)
    x, s = small_roots(f.change_ring(R), bounds, m=3, d=5)[0]
    print(x, s)
    ed = f(x, s)
    assert power_mod(87, ed, N) == 87