Coppersmith’s method is a familiar tool for finding small roots of polynomial congruences.Although this cryptanalysis tool is not currently possible for any cryptography scheme , this topic is worth learning or we can apply for ctf challenge . In RSA challenges, this often means recovering an unknown part of a message or a secret key from partial information. When the problem involves two unknowns, we need to model their relationship with a bivariate polynomial and construct a lattice that captures the bounds on both variables.

In this post, I walk through that construction step by step: choosing shifted polynomials, scaling their coefficients by the root bounds, applying LLL, and extracting candidate roots from the resulting polynomial relations. I also share my SageMath implementation and a small worked example to make the matrix construction easier to follow.

The second half applies this approach to three RSA settings: a stereotyped message with two unknown parts, a small private exponent in the Boneh–Durfee attack, and partial exposure of the private key.

BIVARIATE COPPERSMITH

1. Overview

1.1. Howgrave-Graham

• Problem: Find a small solution $\vert{}x_0\vert{} \le X, \vert{}y_0\vert{} \le Y$ such that $f(x, y) \equiv 0 \pmod M$.

• Theorem: Let $f(x, y) \in \mathbb{Z}[x, y]$ be an irreducible polynomial of degree $\delta$. The target solution $(x_0, y_0)$ is bounded by $\vert{}x_0\vert{} \le X$ and $\vert{}y_0\vert{} \le Y$. Set $W = \max_{i,j} \vert{}f_{i,j}\vert{} X^i Y^j$. If the condition holds: $XY < W^{1/\delta}$ then all solutions $(x_0, y_0)$ can be found in polynomial time in $(\log W, 2^\delta)$.

• Practical root bound: Not dependent on $W^{1/\delta}$, but derived from the inequality: $\det(\mathcal{L}) < M^{s_M} \iff X^{s_x} Y^{s_y} < M^{s_M}$

• Mechanism: Construct a lattice from shifted polynomials containing powers modulo $M^{m-k}$. The goal is to find at least two short linearly independent vectors satisfying the Howgrave-Graham condition: $\Vert{}\mathbf{b}\Vert{} < \frac{M^m}{\sqrt{n}}$ in order to transform them into two algebraic equations over the integers: $h_1(x, y) = 0$ and $h_2(x, y) = 0$.

  • Then use the resultant (or Gröbner basis) to eliminate one variable: $r(x) = \text{resultant}(h_1(x, y), \; h_2(x, y), \; y)$
  • However, LLL only guarantees that the two found vectors are linearly independent, not necessarily algebraically independent. In practice, for RSA-type cryptographic experiments, the algorithm often succeeds with very high probability.

2. Lattice construction method

2.1. Matrix construction rules

• 1-to-1 monomial matching: The matrix must be lower triangular (or full rank). Each newly generated monomial in the columns must correspond to exactly one shifted polynomial whose leading monomial matches the row.
• Replace the variables $x \to Xx$ and $y \to Yy$.

2.2. Parameters $(m, t)$ and matrix size

• Parameter $m$ (modulo power):
  • $m$ is the maximum exponent of $f(x, y)$; larger $m$ implies more monomials in $f(x, y)^m$, so the matrix grows.
  • Larger $m$ yields a wider theoretical bound, but the matrix size grows roughly as $\mathcal{O}(m^2)$.
• Parameter $t$ (extra shifts):
  • The monomials are not only those in $f(x, y)^m$, but also newly generated monomials from $x^i y^j f(x, y)^m$ with $(i+j \le t)$.
  • $t$ is the bound on the shift monomial exponents.
• The matrix dimension $n$ is the total number of unique monomials appearing after collecting all shifted polynomials.

2.3. General formula for shifted polynomials

The family of shifted polynomials modulo $M$ is: $g(x, y) = x^i y^j \cdot [f(x, y)]^k \cdot M^{m-k}$

• $k \in \{0, 1, \dots, m\}$: exponent of the target polynomial.
• $M^{m-k}$: complementary power of the modulus so that at the root $(x_0, y_0)$, every shifted polynomial is divisible by $M^m$.
• $x^i y^j$: shift factor used to generate new leading monomials and fill the matrix.

2.4. Optimization goal

• Objective: To make the LLL algorithm produce vectors satisfying the Howgrave-Graham theorem, we need to reduce the average factor $\det(L)^{1/n}$ as much as possible by tuning the shift parameters $m$ and $t$.

• General optimal inequality: $\det(\mathcal{L}) = M^{S_M} \cdot X^{S_x} \cdot Y^{S_y}$ where $S_M, S_x, S_y$ are the sums of the exponents of $M, X, Y$ on the main diagonal. We need the following inequality to hold for the lattice to work: $\det(\mathcal{L}) < M^{m \cdot n}$ $S_M + S_x \log_M(X) + S_y \log_M(Y) < m \cdot n$

• Optimization idea:

  • We tune $m$, $t$, and the shift pattern so that the left-hand side $S_M + S_x \log_M(X) + S_y \log_M(Y)$ grows more slowly than the right-hand side $m \cdot n$.
  • Possible approaches:
    • Increase $m$: The inequality approaches the theoretical limit as $m \to \infty$, but this may make the matrix much larger.
    • Adjust $t$ and the shift pattern: The parameter $t$ controls the extra shifts. Increasing $t$ may enlarge $n$, but it can also increase $S_M, S_x, S_y$. We can tune $t$ and focus on shifts in the variable with the smaller bound.

2.5 Implement

I implemented bivariate coppersmith via this code

import itertools 


def solve_sys_2_equations_non_trivial(f1,f2): 
    x, y = f1.parent().gens()

    res_x = f1.resultant(f2,y)               
    if not res_x:
        return None 

    roots_x = res_x.univariate_polynomial().roots()

    roots = []
    for rootx in roots_x:
        x_0 = rootx[0]
        if x_0 != 0 :
            PR_y.<y_var> = PolynomialRing(ZZ)
            f1_eval = f1(x_0, y_var)
            y_roots = f1_eval.roots()
            
            for rooty in y_roots:
                y_0 = rooty[0]
                if y_0 != 0 :
                    roots.append((x_0, y_0))
    return roots 




def gen_shift (f,m,target_monos, M):
    shifts = []

    lead_monos = {i:(f^i).lm() for i in range(m+1)}

    for mono in target_monos:
        for i in range(m,-1 ,-1):
            a = mono.degree(x)
            b = mono.degree(y)
            u = lead_monos[i].degree(x)
            v =  lead_monos[i].degree(y)
            if a >=u and b >= v:
                shift_poly = x^(a-u) * y^(b-v) * (f^i) * M ^(m-i)
                shifts.append(shift_poly)
                break 

    return shifts 




def small_root_bivariate(f ,M,bounds , m , extra_shift):
    x, y = f.parent().gens()
    X,Y = bounds[0], bounds[1]


    lead_coe = f.lc()
    if lead_coe != 1:
        try:
            lead_coe_inv = pow(int(lead_coe), -1, int(M))
            f_monic = 0
            for coeff, mono in f:
                new_coeff = (coeff * lead_coe_inv) % M
                if new_coeff > M // 2:
                    new_coeff = new_coeff - M
                f_monic += new_coeff * mono
            f = f_monic
        except:
            print("polymomial can't be monic, should minimize lead coe if it is big")

            
    base_monomials = (f^m).monomials()
    base_shifts = gen_shift(f,m,base_monomials,M)
    

    extra_monomials = []

    for extra_mono in extra_shift:
        extra_poly = f^m * extra_mono
        extra_poly_monos = extra_poly.monomials()
        for mono in extra_poly_monos:
            if mono not in base_monomials and mono not in extra_monomials:
                extra_monomials.append(mono)

    extra_shifts = gen_shift(f,m,extra_monomials,M)
   

    all_monos = base_monomials + extra_monomials
    all_monos = sorted(all_monos, key=lambda m: m.degree())
    all_shifts = base_shifts + extra_shifts
    dim1 = len(all_shifts)
    dim2 = len(all_monos)
    print("Size of lattice: ",dim1 ,"x", dim2)


    M = matrix(ZZ,dim1,dim2)
    
    for r,poly in enumerate(all_shifts):
        for c , mono in enumerate(all_monos):
            i = mono.degree(x)
            j = mono.degree(y)
            coe = poly.monomial_coefficient(mono)
            M[r,c] = coe * X^i * Y^j

    
    M = M.LLL()

    reduced_polys = []
    for row in M:
        poly = 0 
        for c , mono in enumerate(all_monos):
            i = mono.degree(x)
            j = mono.degree(y)
            poly += ZZ(row[c] // (X^i*Y^j)) * mono
        reduced_polys.append(poly)


    for h1, h2 in itertools.combinations(reduced_polys[:4], 2):
        if h1.is_constant() or h2.is_constant():
            continue

        roots = solve_sys_2_equations_non_trivial(h1, h2)
        if roots:
            print('Root exist\n')
            return roots
    
    print('No root found\n')
    return None

3. Example 1: Basic bivariate case ($m = 2, t = 0$)

Consider the polynomial: $f(x, y) = xy + ax + b \equiv 0 \pmod N \quad (\vert{}x_0\vert{} < X, \; \vert{}y_0\vert{} < Y)$ Since $t = 0$, we only consider the monomials that naturally appear in the expansion of $f(x, y)^2$.

Expanding $f^2$ yields exactly 6 monomials: $\mathcal{M} = [1, \quad x, \quad x^2, \quad xy, \quad x^2y, \quad x^2y^2]$

3.2. System of 6 shifted polynomials (rows)

• Layer $k = 0$ ($N^2$):
  • $g_1 = x^0 f^0 N^2 = \mathbf{N^2}$
  • $g_2 = x^1 f^0 N^2 = \mathbf{N^2 x}$
  • $g_3 = x^2 f^0 N^2 = \mathbf{N^2 x^2}$
• Layer $k = 1$ ($N \cdot f$):
  • $g_4 = x^0 f^1 N = N(xy + ax + b) = \mathbf{N xy} + aN x + bN$
  • $g_5 = x^1 f^1 N = Nx(xy + ax + b) = \mathbf{N x^2y} + aN x^2 + bN x$
• Layer $k = 2$ ($f^2$):
  • $g_6 = x^0 f^2 N^0 = (xy + ax + b)^2 = \mathbf{x^2y^2} + 2a x^2y + a^2 x^2 + 2b xy + 2ab x + b^2$

3.3. Lattice matrix $6 \times 6$ (after scaling $Xx, Yy$)

$L = \begin{bmatrix} N^2 & 0 & 0 & 0 & 0 & 0 \\ 0 & X N^2 & 0 & 0 & 0 & 0 \\ 0 & 0 & X^2 N^2 & 0 & 0 & 0 \\ b N & a X N & 0 & X Y N & 0 & 0 \\ 0 & b X N & a X^2 N & 0 & X^2 Y N & 0 \\ b^2 & 2ab X & a^2 X^2 & 2b X Y & 2a X^2 Y & X^2 Y^2 \end{bmatrix}$

• Main diagonal: $[N^2, \; X N^2, \; X^2 N^2, \; X Y N, \; X^2 Y N, \; X^2 Y^2]$
• Determinant: $\det(L) = N^8 \cdot X^8 \cdot Y^4$

RSA stereotyped multipart

Let consider some example problem that we can use bivariate coppersmith . We knew about RSA stereotyped attack , but if there is not only 1 unknonwn part in the message , we need to use multivariate coppersmith

Suppose a RSA system with

• modulus N
• small e
• leak a lot content of the message , and the unknowm parts of the message are small (compared to modulus N in a threshold)

=> We can apply multivariate coppersmith to recover the full leaked message

Example

Let leaked message

m = “my four letter username is XXXX and my secret four digit pin code is YYYY′′

• where the XXXX and YYYY are unknown parts of the message we to recover. We can rewrite m as $m = m′ + 2^{t_x} x + 2^{t_y}y$
  • where x and y represent the unknown parts of the message and m′ represents the known part:
  • The $2^{t_x}$ and $2^{t_y}$ factors exist to capture position of the unknown parts in the message. In this example, $t_x = 42 × 8 = 336$ and since ty starts at the least significant bit, $t_y = 0$.

m′ =“my four letter username is \x00\x00\x00\x00 and my secret four digit pin code is \x00\x00\x00\x00′′

To recover $x$ and $y$ , we construct the modular polynomial $f (x, y)$ of degree e in x and y:

• Then use Coppersmith’s method for multivariate polynomials to recover the small root (x , y ).

I have construct this attack with very small x and y

#from Crypto.Util.number import bytes_to_long
#n = 1455765726967875294222042036176582183046639449319799865257725048632997208206412725060464539896108386151699035861012094022710149721200601647493419532977859797179257773230596249324811
#e = 2
#msg = 'my four letter username is 3456 and my secret four digit pin code is 1105'
#ct = pow(bytes_to_long(msg.encode()), e, n)
#print(ct)


from sage.all import *
load('/home/nhat-anh/crypto/lattice/coppersmith/bivariate_coppersmith.sage')

e = 2
n = 1455765726967875294222042036176582183046639449319799865257725048632997208206412725060464539896108386151699035861012094022710149721200601647493419532977859797179257773230596249324811
ct = 213510694818628662951848929802201599490601145171353386858417186415677154983136586805325904157950263706190078214327641151833402996395736085170015367777503163712669748251914362251734

msg = b'my four letter username is XXXX and my secret four digit pin code is YYYY'

pos_x = len(' and my secret four digit pin code is YYYY')*8
m_leak = int.from_bytes(msg,'big') - int.from_bytes(b'XXXX','big')*2^pos_x - int.from_bytes(b'YYYY','big')


P.<x,y> = PolynomialRing(ZZ)
f = (m_leak + x*2^pos_x + y)^e - ct

bounds = (2^32, 2^32)

roots = small_root_bivariate(f,n,bounds ,3,[y,x])

if roots:
    print("Root :",roots)
    x0,y0 = roots[0][0] , roots[0][1]
    full_msg = int(m_leak + x0*2^pos_x + y0) 
    print(int.to_bytes(full_msg,full_msg.bit_length()//8 +1 ,'big'))

Size of lattice:  36 x 36
Root exist

Root : [(859059510, 825307189)]
b'my four letter username is 3456 and my secret four digit pin code is 1105'

Boneh-Durfee Attack

Wiener (1990): Using the continued fraction technique, Wiener proved that an RSA cryptosystem can be completely broken when the decryption exponent is small: $d < N^{0.25}$

Boneh & Durfee (1999): Using a new approach based on lattice reduction techniques (Lattice Reduction / Coppersmith’s method), Boneh and Durfee extended the attack bound to: $d < N^{0.292}$

Attack

Problem:

• $N = pq$ is an RSA modulus with balanced factors ($p, q$ have the same bit length, $p \approx q \approx \sqrt{N}$).
• $e$ is the public key, and $d$ is the private key satisfying: $ed \equiv 1 \pmod{\phi(N)} \quad \text{with } d < N^\delta$

Polynomial setup

By the definition of congruence, there exists an integer $k$ such that: $ed = 1 + k \cdot \phi(N)$ Substituting $\phi(N) = (p-1)(q-1) = N - (p+q) + 1$: $ed = 1 + k(N - p - q + 1)$

Since $ed$ is divisible by $e$ ($ed \equiv 0 \pmod e$), reduce both sides modulo $e$: $1 + k(N - p - q + 1) \equiv 0 \pmod e$

Because $p$ and $q$ are odd primes:

• $N = pq$ is odd $\implies N + 1$ is even $\implies \dfrac{N+1}{2} \in \mathbb{Z}$.
• $p + q$ is the sum of two odd numbers $\implies p + q$ is even $\implies \dfrac{p+q}{2} \in \mathbb{Z}$.

Factor out $2$: $1 + 2k \left( \frac{N+1}{2} - \frac{p+q}{2} \right) \equiv 0 \pmod e$

Let the known constant $A = \dfrac{N+1}{2}$. Then the equation becomes a bivariate polynomial root-finding problem: $f(x, y) = 1 + 2x(A - y) \equiv 0 \pmod e$

The integer solution we need is: $(x_0, y_0) = \left( k, \; \frac{p+q}{2} \right)$

Bounds

$k = \frac{ed - 1}{\phi(N)} < \frac{ed}{\phi(N)}$ For the standard case $e < \phi(N)$, we have $k < d$. Therefore: $|x_0| = k < X = N^\delta$

Since the two prime factors are balanced ($p \approx q \approx \sqrt{N}$), we have: $|y_0| = \frac{p+q}{2} < Y \approx N^{0.5}$

LLL and solving the system

• Find the solution using lattice methods: apply the two-variable Coppersmith method with shifted polynomial systems ($x$-shifts and $y$-shifts) modulo $e^m$.

  • In Boneh-Durfee’s work, it is shown that using $y$-shifts ($y, \ldots, y^t$) can enlarge the root bound.
  • The same study also gives the optimal value of $t$: $t \approx (1 - 2\delta)m$
  • One can run the Coppersmith algorithm with increasing modulus $m$ until the root is found; the algorithm guarantees finding the integer solution $(x_0, y_0) = (k, s)$ in polynomial time if $\delta < 0.292$

• Recover the key $d$ after finding $(k, s)$ There are two equivalent ways to recover the private key:

  $ed = 1 + 2k(A - s) \implies d = \frac{1 + 2k\left(\frac{N+1}{2} - s\right)}{e}$

Example

I have create a toy ctf challenge about the Boneh-Durfee attack

• The challenge:

  
  from Crypto.Util.number import getPrime, inverse, bytes_to_long
  import random
  
  while True:
      p = getPrime(256)
      q = getPrime(256)
      if 0.9 < p/q < 1.1:
          break
  n = p*q
  print("n =", n)
  phi = (p-1)*(q-1)
  
  R = RealField(512)
  
  while True: 
      d = random.randint(int(R(n)^0.25), int(R(n)^0.27))
      if gcd(d,phi) ==1: 
          e = inverse(d, phi)
          print("e =", e)
          break 
  
  m = b'LTNA{b1v4r1at3_C0pperSm1th}'
  ct = pow(bytes_to_long(m), e, n)
  print("ct =", ct)

• The solution

  from Crypto.Util.number import long_to_bytes 
  load('bivariate_coppersmith.sage')
  
  n = 7231784656254883250070144994499124584107632716753742081690327324192790224881727213289232394149250983798248480282721659908052976120975159475848896467902137
  e = 1679059439441963415529389383875030860272633267370001412291746837734438272855590503547653168455264858393609768228046928485997033002304953651224626430513779
  ct = 4361287623317308583715975471074068285530921515353045304187828853936145692776409661620168196922237957137883723956335629774878881760908412405505224838532385
  
  RF = RealField(256)
  bounds = [int(RF(n)^0.27), int(RF(n)^0.5)]
  
  A = (n+1)//2 
  PR.<x,y> = PolynomialRing(ZZ, 2)
  f = -2*x*y + 2*A*x + 1
  
  
  delta = 0.27
  
  for m in range(2,10):
      t= int(m*(1-2*delta))
      extra_shift = [y^i for i in range(1,t+1)]
  
      roots = small_root_bivariate(f ,e,bounds , m , extra_shift)
  
      if roots:
          x0, y0 = roots[0][0], roots[0][1]
          d = (1+2*x0*((n+1)//2 - y0))// e 
  
          print("d = ",d)
          flag = long_to_bytes(pow(ct,d,n))
          print("flag :",flag)
          break 

  d =  209033656159098247762466303108200858417339
  flag : b'LTNA{b1v4r1at3_C0pperSm1th}'

Partial Key Exposure Attack

Partial key exposure attack’s general approach is quite similar to that of the Boneh-Durfee attack.

Problem

RSA system

• balanced modulus $N = pq$ with ($p \approx q$)

• public exponent $e$, and the corresponding private exponent $d$.

• Furthermore, we are given that $d < N^\delta$ as well as $(\delta - \gamma) \log_2 N$ MSBs of $d$ => we know an integer $d_0$ such that $|d - d_0| < N^\gamma$.

Solve

If $\delta < 0.292$, then we can directly apply the Boneh-Durfee attack without using the additional partial knowledge of $d$.

If $\delta > 0.292$, it is still helpful to revisit the reformulation of the RSA key equation used in the attack:

$1 + 2k \left( \frac{N + 1}{2} - \frac{p + q}{2} \right) \equiv 0 \pmod e$

The key idea is that we can use the known bits of $d$ to approximate $k$ so that the bounds on the roots that we need to solve for are lowered.

• Let $d = d_0 + d_1$ where $d_1$ is the remainder of the unknown bits of $d$.

$k = \frac{ed - 1}{\varphi(N)} = \frac{e(d_0 + d_1) - 1}{\varphi(N)} \approx \frac{ed_0}{N}$

• $k_0 = \left\lfloor \frac{ed_0}{N} \right\rfloor$ is a good approximation for $k$. It is that this approximation satisfies $|k - k_0| < 4N^\lambda$ with $\lambda = \max\left\{\gamma, \delta - \frac{1}{2}\right\}$.
• Therefore, writing $k = k_0 + k_1$ where $k_1$ is the remainder of the unknown bits of $k$, we find that the polynomial

$f(x, y) = 1 + 2(k_0 + x) \left( \frac{N + 1}{2} - y \right) \pmod e$

• $f(x, y)$ has the “small” root $(x, y) = \left(k_1, \frac{p + q}{2}\right)$ for appropriate $\delta$ and $\gamma$.

Example

I have create a toy ctf challenge about the Boneh-Durfee attack

• The challenge:

  from Crypto.Util.number import getPrime , long_to_bytes , bytes_to_long , inverse
  
  
  while True:
      p = getPrime(256)
      q = getPrime(256)
      if 0.9 < p/q < 1.1:
          break
  n = p*q
  print("n =", n)
  phi = (p-1)*(q-1)
  
  
  while True: 
      d = getPrime(320)
      e = inverse(d, phi)
      if e.bit_length() == n.bit_length() and e < n:
          print("e =", e)
          d_leak = d >> 70
          print("d_leak =", d_leak)
          break 
  
  
  flag = b'LTNA{P4rt1Al_k3y_exP0sur3}' 
  ct = pow(bytes_to_long(flag),e,n)
  print("ct =",ct)
  
  
  
  """
  n = 9533706189204254935659299342308408208072430877029789122313277458953249322469389045357064238990459147749307149206576043308797716161128494068245409126585863
  e = 7708858996381546599613326302093374442396107334245800056896760670090587468486003428457987002295831136473104728695953066752332518276482112750045111725973991
  d_leak = 1090150921049954754060630900546022989625881542269502482570883961432553998392
  ct = 2546148024525155120114602408888547250891277919807461327345001225843804052055433979861218949306846252464948048273131967790896726974775561256126448476936679
  """

• The solution

  from Crypto.Util.number import long_to_bytes 
  load("bivariate_coppersmith.sage")
  
  
  n = 9533706189204254935659299342308408208072430877029789122313277458953249322469389045357064238990459147749307149206576043308797716161128494068245409126585863
  e = 7708858996381546599613326302093374442396107334245800056896760670090587468486003428457987002295831136473104728695953066752332518276482112750045111725973991
  d_leak = 1090150921049954754060630900546022989625881542269502482570883961432553998392
  ct = 2546148024525155120114602408888547250891277919807461327345001225843804052055433979861218949306846252464948048273131967790896726974775561256126448476936679
  
  bounds = [2^72, 2^256]
  
  A = (n+1)//2 
  
  d_leak = d_leak << 70 
  k0 = (e * d_leak) // n
  
  PR.<x,y> = PolynomialRing(ZZ, 2)
  f = 1 + 2*(k0 + x)*(A - y)
  
  roots = small_root_bivariate(f, e, bounds, 5, [x, x^2])
  
  if roots: 
      x0, y0 = roots[0][0], roots[0][1]
      k = k0 + x0
      d = (1+2*k*((n+1)//2 - y0))// e 
  
      print("d =",d)
      flag = long_to_bytes(pow(ct,d,n))
      print("flag :",flag)

  Size of lattice:  48 x 48
  Root exist
  
  d = 1287023042708944777222463971824550736964155950203023319856968787964858150943107946377689028308967
  flag : b'LTNA{P4rt1Al_k3y_exP0sur3}'