Continued Fractions And Lattice Sieving

A newsletter about applied math and theoretical computer science. Subscribe!

This is part of our Practical Number Field Sieve for Programmers series:

Part 1: Discrete Logarithms and the Index Calculus Solution.

Part 2: Solving Index Calculus Equations over Integers and Finite Fields.

Part 3: Computation of Discrete Logarithms in Algebraic Number Fields.

Part 4: 2 Dimensional Lattice Basis Reduction.

Part 5: Individual Reduction Phase and Logarithm Collection.

Part 6: Lattice Sieving and Special Q Descent.

Part 7: Continued Fractions and Fast Lattice Sieving.

Code is available on Colab.

1.0 Paper Introduction

Continued Fractions and Lattice Sieving (Franke & Kleinjung, 2025) improves on both Lenstra’s sieve and the Golliver, Lenstra & McCurley (GLM) lattice sieve encountered in Part 6.

(Franke & Kleinjung, 2025) address the core challenge in the GLM sieve: finding the next valid lattice point is somewhat tedious.

The authors depend on continued fraction coefficients to yield the best rational approximation to a real number (Krishnan, 2016)

This construction permits easy jumps among valid lattice points as summarized in (Zimmerman, 2009)

1.1 Continued Fractions and Lattice Point Enumeration

As we saw earlier, a special-q pair is a pair (q,ρ) such that q is a prime and ρ is the root of the number field sieve polynomial used in our finite field.

We start from the special-q basis constructed in section 1.2 of Part 6:

We want to find integers a, b such that:

Furthermore, we want our basis vectors to satisfy the three conditions:

(Franke & Kleinjung, 2025) use the fact that the continued fraction expansion of r/z gives the best (pk, qk) pairs such that:

The next section finds r/z as demonstrated in (Franke & Kleinjung, 2025).

1.2 Algorithm 1: Generating Franke-Kleinjung Basis

Code is available on Colab.

First, we define a basis, and use the algorithm below to reduce it:

Python code resembles

def Algorithm1(I, u0, u1):
  u1 = [u1[0] - u0[0], u1[1] - u0[1]]
  while abs(u1[0]) >= I:
    a = math.floor(-u0[0] / u1[0])
    u0 = [u0[0] + a*u1[0], u0[1] + a*u1[1]]
    u0, u1 = u1, u0
    a = math.floor((abs(u0[0]) - I) / abs(u1[0])) + 1
    u0 = [u0[0] + a*u1[0], u0[1] + a*u1[1]]
  return tuple(u0), tuple(u1)

For the observant, Algorithm 1 involves finding continued fraction coefficients, hence the paper’s name.

1.3 Algorithm 2: Enumerating All Lattice Points

The basis constructed in Algorithm 1 has the special property that the second basis vector increases monotonically.

This permits us start from the zero vector (0,0) and walk the lattice by recursively adding u0 and u1.

In Python, this resembles:

def NEXTFK2(I, u0, u1, p):
  p0, p1 = p
  half = I // 2
  if -half <= p0 + u0[0]:
      return (p0 + u0[0], p1 + u0[1])
  if p0 + u1[0] < half:
      return (p0 + u1[0], p1 + u1[1])
  return (p0 + u0[0] + u1[0], p1 + u0[1] + u1[1])

def EnumeratePoints(I, u0, u1, maxPoints=10):
  """Enumerate all points in L ∩ H starting from (0,0)"""
  p = (0, 0)
  points = [p]
  for i in range(maxPoints - 1):
      p = NEXTFK2(I, u0, u1, p)
      points.append(p)
  return points

2.0 Implementation Notes

This section ends with a concrete demonstration of the sieve.

2.1.1 Other Class Numbers

We use this table to find norms and roots for other Heegner numbers:

For instance, for 43, one evaluates the ring modulo p to find the root:

2.1.2 Algebraic Sieve

Our number field sieve builds on this esoteric fact from Part 5:

That is, we use rational reconstruction from Section 3 to find T and V such that:

We follow Part 3 and work with a 191-bit prime alongside two special-q’s that split over the Heegner number -67:

p = 2606809712135760185843507056585091170888781801840177208389
q0 = 2941869594359
q1 = 17
negativeUnderRoot = -67

We use the chinese remainder theorem to find z where q=q0 and r=q1 and rho is the number field root:

Next, we simply walk to enumerate points. The exact code resembles:

from sympy.ntheory.residue_ntheory import sqrt_mod
from sympy.ntheory.modular import crt
from sympy.ntheory import factorint
import math

def dot(u, v):
  return u[0]*v[0] + u[1]*v[1]

def norm2(v):
  return dot(v, v)

def GaussReduce(b1, b2):
  while True:
    if norm2(b2) < norm2(b1):
      b1, b2 = b2, b1
    m = round(dot(b1, b2) / norm2(b1))
    if m == 0:
      return b1, b2

    b2 = (b2[0] - m * b1[0],b2[1] - m * b1[1])

def Algorithm1(I, u0, u1):
  u1 = [u1[0] - u0[0], u1[1] - u0[1]]
  while abs(u1[0]) >= I:
    a = math.floor(-u0[0] / u1[0])
    u0 = [u0[0] + a*u1[0], u0[1] + a*u1[1]]
    u0, u1 = u1, u0
    a = math.floor((abs(u0[0]) - I) / abs(u1[0])) + 1
    u0 = [u0[0] + a*u1[0], u0[1] + a*u1[1]]
  return tuple(u0), tuple(u1)

def FindNorm(a, b, negativeUnderRoot):
  if negativeUnderRoot % 4 == 1:
    #alpha = (1 + sqrt(negativeUnderRoot)) / 2
    D14 = (negativeUnderRoot - 1) // 4
    return a*a + a*b - D14*b*b
  else:
    #alpha = sqrt(negativeUnderRoot)
    return a*a - negativeUnderRoot*b*b

def FindRoot(q, negativeUnderRoot):
  s = sqrt_mod(negativeUnderRoot, q)
  assert s is not None
  if negativeUnderRoot % 4 == 1:
    #print(f"Gere:{negativeUnderRoot}, {s}, q:{q}")
    #alpha = (1 + sqrt(negativeUnderRoot)) / 2
    return ((s + 1) * pow(2, -1, q)) % q
  else:
    #alpha = sqrt(negativeUnderRoot)
    return s % q

def NEXTFK2(I, u0, u1, p):
  p0, p1 = p
  half = I // 2
  if -half <= p0 + u0[0]:
      return (p0 + u0[0], p1 + u0[1])
  if p0 + u1[0] < half:
      return (p0 + u1[0], p1 + u1[1])
  return (p0 + u0[0] + u1[0], p1 + u0[1] + u1[1])

def EnumeratePoints(I, u0, u1, maxPoints=10):
  """Enumerate all points in L ∩ H starting from (0,0)"""
  p = (0, 0)
  points = [p]
  for i in range(maxPoints - 1):
      p = NEXTFK2(I, u0, u1, p)
      points.append(p)
  return points

p = 2606809712135760185843507056585091170888781801840177208389
q0 = 2941869594359
q1 = 17
negativeUnderRoot = -67
maxPoints=2
I = 2 ** 6
rho_q0 = FindRoot(q0,negativeUnderRoot)
rho_q1 = FindRoot(q1,negativeUnderRoot)

rho_p = FindRoot(p,negativeUnderRoot)
b1 = (p, 0)
b2 = (-rho_p, 1)
r1, r2 = GaussReduce(b1, b2)
T = r1[0]
V = r1[1]
test = (T + V * rho_p) % p
print(f"test TV: {test}, T: {T}, V:{V}")
print(r1)
print(r2)
z, r = crt([q0, q1], [-rho_q0, -rho_q1])
z = z % r
u0 = (r, 0)
u1 = (z, 1)
result = Algorithm1(I, u0, u1)
points = EnumeratePoints(I, result[1], result[0], maxPoints)
print(f"\nStart basis: u0={u0}, u1={u1}")
print(f"rho_q0 = {rho_q0}\nrho_q1 = {rho_q1}\nrho_p = {rho_p}\nr = {r}\nz = {z}")
print(f"\nFinal basis: u0={result[0]}, u1={result[1]}")
for i, (a0, a1) in enumerate(points):
    if a0 == 0 and a1 == 0:
        continue
    N = FindNorm(a0, a1, negativeUnderRoot)
    factors = factorint(N)
    factorString = " × ".join([f"{k}^{v}" if v > 1 else str(k) for k, v in factors.items()])
    lhs = (a0*V - a1 * T) % p
    rhs = V*(a0 + a1 * rho_p) % p
    q0Valid = (N % q0 == 0)
    q1Valid = (N % q1 == 0)
    print(f"{i:6d}: ({a0}, {a1}),{q0Valid and q1Valid}: N = {N}")
    print(f"     = {factorString}")
    if(lhs != 0):
      lhsSize = math.log(abs(a0*V - a1 * T), 2)
      print(f"cV-dT map: {lhs == rhs} : lhsSize: {lhsSize:.2f} : \n")

2.1.3 Rational Sieve

Section 2.1.2 demonstrated lattice enumeration on the algebraic side. This section demonstrates the walk on the rational side.

We desire to find c, d pairs such that the LHS and RHS are smooth in:

This is achieved by sieving over x, y multiples of the original basis. That is, we desire to find x and y such that:

Furthermore, c and d satisfy the relation below modulo the small primes mi:

We follow the sieve approach from Part 5, to construct a rational sieve where we fix c and find a good d. This construction permits us identify great x, y lattice points that are smooth on both the rational and algebraic side.

def RationalSieve(basis0, basis1, T, V, primes, X, Y):
  candidates = []
  primes=[2,3,5,7,11,13,17,19,23,29,53,257]
  for prime in primes:
    det = (u0[0] * u1[1] - u1[0] * u0[1]) % prime
    if(T % prime == 0):
      continue
    if(det % prime == 0):
      continue
    TInverse = pow(T, -1, prime)
    detInverse = pow(det, -1, prime)
    k = (V * TInverse) % prime
    xStep = ((u1[1] - k*u1[0]) * detInverse) % prime
    yStep = ((k*u0[0] - u0[1]) * detInverse) % prime
    slope = (yStep * pow(xStep, -1, prime)) % prime
    #slope says invcrease in x by 1, y increases by slope
    for c in range(1, prime):
      #d=cVTInverse mod prime
      d = (c * V * pow(T, -1, prime)) % prime
      x = ((c * u1[1] - d * u1[0]) * pow(det, -1, prime)) % prime
      y = ((d * u0[0] - c * u0[1]) * pow(det, -1, prime)) % prime
      print(f"det:{det:3d}, x:{x:3d} , y:{y:3d}, slope: {slope} mod {prime}")
    print("")
  return candidates
u0 = result[1]
u1 = result[0]
print(f"\nFinal basis: u0={u0}, u1={u1}")
candidates = RationalSieve(u0, u1,T, V,primes=[2,257],X=100,Y=100)

for c, d in candidates:
  R = c*V - d*T
  print(c, d, R)

Subscribe if you made it this far!

References

Franke, J & Kleinjung T. (2025). Continued Fractions and Lattice Sieving. Hyperelliptic: Proceedings of SHARCS 2005. PDF Link.

Krishnan, G. (2016).Continued Fractions. Cornell University. PDF.

Zimmerman, P. (2009). CADO-NFS: An Implementation of The Number Field Sieve. Loria INRIA. PDF Link.