Bivariate Coppersmith Algorithm

edit

Given the congruences:

(a*n-a+2) mod 83 = 0
(a*n-1884*a+119) mod 1683 = 0
(a*n-731*a+649) mod 5923 = 0
(a*n-6654*a+6572) mod 53307 = 0
(a*n-46652*a+46570) mod 373291 = 0
(a*n-46653*a+46571) mod 373299 = 0

and

|n|<|a|<27

ChatGPT wrote this for me: (Please excuse me for using AI)

# SageMath
# Reticolo bivariato LLL per le sei congruenze
# an + c_i*a + d_i == 0 (mod N_i)
#
# Bound: |n| < |a| < 27
#
# Nota: e' un metodo di tipo Coppersmith/Howgrave-Graham
# euristico per il sistema combinato, non una garanzia
# teorica di recupero delle radici.

from sage.all import *

# --------------------------------------------------
# 1. Dati
# --------------------------------------------------

Ns = [
    ZZ(83),
    ZZ(1683),
    ZZ(5923),
    ZZ(53307),
    ZZ(373291),
    ZZ(373299)
]

cs = [
    ZZ(-1),
    ZZ(-1884),
    ZZ(-731),
    ZZ(-6654),
    ZZ(-46652),
    ZZ(-46653)
]

ds = [
    ZZ(2),
    ZZ(119),
    ZZ(649),
    ZZ(6572),
    ZZ(46570),
    ZZ(46571)
]

X = ZZ(26)
Y = ZZ(25)

# Parametri del reticolo: aumentabili
t = 3

# --------------------------------------------------
# 2. Polinomi bivariati
# --------------------------------------------------

R.<a,n> = PolynomialRing(ZZ, order='lex')

f = [
    a*n + cs[i]*a + ds[i]
    for i in range(len(Ns))
]

print("Polinomi originali:")
for i in range(len(Ns)):
    print("f_%d =" % i, f[i])
    print("modulo", Ns[i])

# --------------------------------------------------
# 3. Combinazione modulare
# --------------------------------------------------
# P e' il prodotto dei moduli.
# Se tutte le congruenze sono soddisfatte, allora
# G(a,n) == 0 mod P.
#
# G = sum (P/N_i)*f_i.
#
# La condizione inversa non e' garantita:
# ogni candidato sara' verificato su tutte le f_i.

P = prod(Ns)

G = sum(
    (P // Ns[i]) * f[i]
    for i in range(len(Ns))
)

print("\nP =", P)
print("Grado totale di G =", G.total_degree())

# --------------------------------------------------
# 4. Generazione degli shift
# --------------------------------------------------
# Si usano shift monomiali di G e shift multipli di P.
# La base viene scalata con X^i * Y^j per incorporare
# i bound delle radici.

def exp_pairs(D):
    return [
        (i, j)
        for total in range(D + 1)
        for i in range(total + 1)
        for j in [total - i]
    ]

shifts = []

for i, j in exp_pairs(t):
    shifts.append(a**i * n**j * G)

for i, j in exp_pairs(t + 2):
    shifts.append(P * a**i * n**j)

# Grado massimo dei polinomi di base
D = t + 4
mons_exp = exp_pairs(D)

# Vettori dei coefficienti scalati
rows = []

for h in shifts:
    row = []
    for i, j in mons_exp:
        coeff = h.monomial_coefficient(a**i * n**j)
        row.append(ZZ(coeff) * X**i * Y**j)
    if any(v != 0 for v in row):
        rows.append(row)

M = matrix(ZZ, rows)

print("\nDimensioni della matrice iniziale:", M.nrows(), "x", M.ncols())

# --------------------------------------------------
# 5. Riduzione LLL
# --------------------------------------------------

L = M.LLL()

print("Dimensioni della base ridotta:", L.nrows(), "x", L.ncols())

# Ricostruisce i polinomi dai vettori ridotti,
# rimuovendo la scalatura X^i*Y^j.

S.<u,v> = PolynomialRing(QQ, order='lex')

reduced_polys = []

for row in L.rows():
    h = S.zero()

    for k, (i, j) in enumerate(mons_exp):
        scale = X**i * Y**j
        if scale != 0 and row[k] != 0:
            h += QQ(row[k]) / scale * u**i * v**j

    if h != 0 and h not in reduced_polys:
        reduced_polys.append(h)

print("\nPolinomi non nulli ottenuti:", len(reduced_polys))

for i, h in enumerate(reduced_polys[:8]):
    print("\nh_%d =" % i, h)

# --------------------------------------------------
# 6. Risultanti per ricavare candidati per a
# --------------------------------------------------

T.<z> = PolynomialRing(QQ)

candidate_a = set()

# Cerca coppie di polinomi con risultante non nullo.
for i in range(len(reduced_polys)):
    for j in range(i + 1, len(reduced_polys)):
        h1 = reduced_polys[i]
        h2 = reduced_polys[j]

        try:
            res = h1.resultant(h2, v)
        except Exception:
            continue

        if res == 0:
            continue

        # Il risultante deve dipendere da u.
        try:
            rz = T(res(u=z))
        except Exception:
            continue

        if rz == 0 or rz.degree() <= 0:
            continue

        # Le radici intere sono candidate per a.
        try:
            for root, mult in rz.roots(ZZ):
                if abs(root) < 27 and root != 0:
                    candidate_a.add(ZZ(root))
        except Exception:
            pass

print("\nCandidati a dai risultanti:", sorted(candidate_a))

# --------------------------------------------------
# 7. Verifica rigorosa nelle congruenze originali
# --------------------------------------------------

def verifica(aa, nn):
    if aa == 0:
        return False
    if not (abs(nn) < abs(aa) < 27):
        return False

    return all(
        (aa*nn + cs[i]*aa + ds[i]) % Ns[i] == 0
        for i in range(len(Ns))
    )

soluzioni = set()

# Per ogni candidato a, n e' cercato solo nel bound
# ristretto. Questa fase filtra i candidati LLL.
for aa in candidate_a:
    for nn in range(-25, 26):
        if verifica(aa, ZZ(nn)):
            soluzioni.add((aa, ZZ(nn)))

print("\nSoluzioni certificate trovate dal reticolo:")
for aa, nn in sorted(soluzioni):
    print("a =", aa, ", n =", nn)

# --------------------------------------------------
# 8. Controllo indipendente entro i bound
# --------------------------------------------------
# Utile per distinguere un fallimento del reticolo
# da una vera assenza di soluzioni.

print("\nControllo esaustivo indipendente:")
soluzioni_esatte = []

for aa in range(-26, 27):
    if aa == 0:
        continue
    for nn in range(-25, 26):
        if abs(nn) < abs(aa) and verifica(ZZ(aa), ZZ(nn)):
            soluzioni_esatte.append((ZZ(aa), ZZ(nn)))

print(soluzioni_esatte)

OUTPUT:

Polinomi originali:
f_0 = a*n - a + 2
modulo 83
f_1 = a*n - 1884*a + 119
modulo 1683
f_2 = a*n - 731*a + 649
modulo 5923
f_3 = a*n - 6654*a + 6572
modulo 53307
f_4 = a*n - 46652*a + 46570
modulo 373291
f_5 = a*n - 46653*a + 46571
modulo 373299

P = 6145999617209997401439561
Grado totale di G = 2

Dimensioni della matrice iniziale: 31 x 36
Dimensioni della base ridotta: 31 x 36

Polinomi non nulli ottenuti: 21

h_0 = -1279368*u*v - 11994075*u + 116902251

h_1 = 1279368*u*v^2 + 1759131*u*v - 95952600*u - 116902251*v + 935218008

h_2 = -1279368*u^2*v - 11994075*u^2 - 1279368*u*v + 104908176*u + 116902251

h_3 = 1279368*u*v^3 + 479763*u*v^2 - 14552811*u*v + 875567475*u - 116902251*v^2 + 1052120259*v - 8533864323

h_4 = 1279368*u^2*v^2 + 1759131*u^2*v - 95952600*u^2 - 28625859*u*v + 1762809183*u - 8066255319

h_5 = -1279368*u^3*v - 11994075*u^3 - 1279368*u^2*v + 104908176*u^2 - 1279368*u*v + 104908176*u + 116902251

h_6 = -1279368*u^2*v - 11994075*u^2 + 1279368*u*v^4 + 479763*u*v^3 - 4317867*u*v^2 + 132254667*u*v - 7751210949*u - 116902251*v^3 + 1052120259*v^2 - 9469082331*v + 76687876656

h_7 = -1279368*u^2*v^3 - 479763*u^2*v^2 + 18390915*u^2*v - 839585250*u^2 + 8155971*u*v^2 - 391806450*u*v + 23931378045*u + 9936691335*v - 153492655563

Candidati a dai risultanti: [-1, 17]

Soluzioni certificate trovate dal reticolo:
a = 17 , n = -4

Controllo esaustivo indipendente:
[(17, -4)]

This is exactly what I was looking for. Is this implementation correct? How could it be improved? What is the computational cost?