Simultaneous Diophantine Approximation with a Lattice

Everyone knows that $355/113$ is a remarkably good approximation of $\pi$ . The simultaneous version of the question asks for a single denominator that works for several numbers at once: an integer $q$ such that $q\alpha_1, \dots, q\alpha_n$ are all close to integers. Continued fractions answer the one-number case, but there is no continued fraction for several numbers. A lattice, built the right way, answers both. This guide builds that lattice, chooses its one free weight, and runs it from one number up to sixty with g6k-rs.

The problem and Dirichlet’s benchmark

Given reals $\alpha_1, \dots, \alpha_n$ and a bound $Q$ , find integers $1 \le q \le Q$ and $p_1, \dots, p_n$ that make every error $|q\alpha_i - p_i|$ small. How small can we hope for? Dirichlet’s theorem, a pigeonhole argument, gives the benchmark:

Dirichlet’s simultaneous approximation theorem

For any reals $\alpha_1, \dots, \alpha_n$ and any integer $Q \ge 1$ there is an integer $1 \le q \le Q$ with

\max_i \, |q\alpha_i - p_i| \;\le\; Q^{-1/n} \quad\text{for suitable integers } p_i.

Equivalently, the quality $\max_i |q\alpha_i - p_i| \cdot q^{1/n}$ can be made at most 1. We use this number throughout: below 1 means “as good as Dirichlet promises”.

With one number the error falls like $1/Q$ ; with ten numbers only like $Q^{-1/10}$ . Sharing a denominator is expensive, which is why the problem is interesting. Dirichlet’s proof is not constructive: it tells you $q$ exists, not how to find it. Trying every $q$ up to $Q = 10^{30}$ is out of the question. The lattice finds a good one in a fraction of a second.

Building the lattice

Follow the embedding recipe. The small object is the vector of errors $(q\alpha_i - p_i)_i$ , together with the requirement that $q$ is not too large. Each error is a linear expression in the unknown integers $q, p_1, \dots, p_n$ , so it fits in a lattice directly. Real numbers cannot go into an integer basis, so first scale by a large integer $M$ and truncate: $A_i = \lfloor \alpha_i M \rfloor$ . Rows are basis vectors:

B = \begin{pmatrix} W & A_1 & A_2 & \cdots & A_n \\ 0 & M & 0 & \cdots & 0 \\ 0 & 0 & M & \cdots & 0 \\ \vdots & & & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & M \end{pmatrix}

Take $q$ copies of the first row and subtract $p_i$ copies of row $i+1$ :

q\,b_0 - \sum_i p_i\, b_i \;=\; \bigl(\, qW,\; qA_1 - p_1M,\; \dots,\; qA_n - p_nM \,\bigr) \;\approx\; \bigl(\, qW,\; M(q\alpha_1 - p_1),\; \dots \bigr).

The first coordinate measures the size of $q$ ; the others measure the errors, magnified by $M$ . A short lattice vector is therefore a small $q$ with small errors, and you read the answer straight off it: $q = b_0/W$ and $p_i = \operatorname{round}(q\alpha_i)$ . The dimension is $n+1$ , and the basis is triangular, so

\det(L) = W \cdot M^{\,n}.

Choosing the weight W

$M$ only sets the precision. The weight $W$ is the real design choice, because reduction minimizes the Euclidean length and does not know which coordinate you care about. A Dirichlet-quality answer has $q \approx Q$ and errors near $Q^{-1/n}$ , so its coordinates are about $QW$ and $MQ^{-1/n}$ . Make them equal:

QW = M Q^{-1/n} \quad\Longrightarrow\quad W = M \, Q^{-(n+1)/n}.

Now the target is balanced: every coordinate is roughly $MQ^{-1/n}$ , and the one free knob, $Q$ , says how large a denominator you are willing to accept. The two ways to get it wrong:

$W$ too small. Large $q$ costs nothing, so reduction returns a huge denominator with tiny errors, far beyond your bound.
$W$ too large. Any $q \ne 0$ becomes expensive and the shortest vectors drift toward small $q$ with poor errors. At the extreme, $W \ge M$ , the rows $(0,\dots,M,\dots,0)$ win and carry $q = 0$ , which answers nothing.

What the Gaussian heuristic predicts

With the balanced weight, $\det(L)^{1/(n+1)} = MQ^{-1/n}$ , so the Gaussian heuristic estimates the shortest vector at

\mathrm{gh}(L) \approx \sqrt{\tfrac{n+1}{2\pi e}}\; M Q^{-1/n}.

A Dirichlet-quality vector has every coordinate at most $MQ^{-1/n}$ , so its length is up to $\sqrt{n+1}\,MQ^{-1/n}$ , which is $\sqrt{2\pi e} \approx 4.13$ times $\mathrm{gh}(L)$ . The planting ratio $\mathrm{gh}/t$ is about 0.24, well below 1. That is the point where this problem differs from most embeddings: there is no single planted answer to isolate. The lattice holds many vectors of about this length, and every one of them is a good simultaneous approximation. You only need a short vector, not the shortest.

The heuristic also predicts the quality you should see. A vector of length $\mathrm{gh}(L)$ has typical coordinates of size $MQ^{-1/n}/\sqrt{2\pi e}$ , which means quality near $1/\sqrt{2\pi e} \approx 0.24$ in low dimension. Reduction does not reach the shortest vector exactly: LLL returns one about $\delta^{\,n+1}$ times longer than $\det(L)^{1/(n+1)}$ , so quality should creep upward with $n$ , and BKZ should slow the creep. Both predictions are heuristics, and the runs below test them. The only obvious competitors, the $q = 0$ vectors, have length at least $M$ , far above $MQ^{-1/n}$ whenever $Q > 1$ .

Precision: why exact integers

Truncating $\alpha_i M$ changes each $A_i$ by less than 1, so each error coordinate is off by less than $q$ . That is harmless when $M$ is large enough that $q \ll MQ^{-1/n}$ ; the code keeps 12 extra digits for this. The consequence is that the entries grow quickly: sixty numbers with $Q = 10^{60}$ need $M = 10^{73}$ , far past the 53-bit mantissa of an f64. So the code uses the exact MpzLatticeBasis of g6k-rs, which keeps every entry as a GMP integer, and it computes $\lfloor\sqrt{p}\cdot 10^k\rfloor$ exactly with an integer square root instead of parsing floating-point values.

The code

The program below is complete and runs as written against g6k-rs = "0.1" (version 0.1.0) plus the rug crate for big integers. The function approximate builds the basis above with the balanced weight, runs LLL (and optionally BKZ), and reads the best $(q, p_1, \dots, p_n)$ off the reduced rows.

// Simultaneous Diophantine approximation with g6k-rs (exact integer arithmetic).
// Cargo.toml:
//   g6k-rs = "0.1"
//   rug = { version = "1.26", default-features = false, features = ["integer", "float"] }
//
// Rows are basis vectors (g6k-rs convention).
//   row 0      = ( W, A_1, A_2, ..., A_n )    with A_i = floor(alpha_i * M)
//   row i (>0) = ( 0, ..., M, ..., 0 )        M in column i
// q*row0 - sum p_i*row_i = ( q*W, q*A_1 - p_1*M, ..., q*A_n - p_n*M ).
use g6k_rs::{MpzLatticeBasis, MpzLllOptions, mpz_bkz, mpz_lll_reduce_with_options};
use rug::Integer;

/// floor(x * 10^k) for a decimal string x with at least k fractional digits.
fn from_decimal(x: &str, k: usize) -> Integer {
    let (int, frac) = x.split_once('.').unwrap();
    format!("{int}{}", &frac[..k]).parse().unwrap()
}

fn pow10(k: usize) -> Integer {
    Integer::from(Integer::u_pow_u(10, k as u32))
}

/// floor(sqrt(p) * 10^k), computed exactly.
fn sqrt_scaled(p: u32, k: usize) -> Integer {
    (Integer::from(p) * pow10(2 * k)).sqrt()
}

struct Approx {
    q: Integer,
    worst: f64, // max_i |q*alpha_i - p_i|
}

/// Dirichlet guarantees some q <= Q with worst < Q^(-1/n), i.e. worst * q^(1/n) < 1.
fn quality(a: &Approx, n: usize) -> f64 {
    a.worst * a.q.to_f64().powf(1.0 / n as f64)
}

/// alphas(k) must return floor(alpha_i * 10^k). Targets q around Q = 10^qe.
/// block = 0: LLL only; block >= 2: LLL then BKZ with that block size.
fn approximate(n: usize, alphas: impl Fn(usize) -> Vec<Integer>, qe: u32, block: usize) -> Approx {
    // Balance q*W against M*|q*alpha - p| ~ M*Q^(-1/n): W = M / Q^((n+1)/n).
    // 12 guard digits keep W a large integer and the truncation error of A_i negligible.
    let e = qe as f64 * (n as f64 + 1.0) / n as f64;
    let k = e.ceil() as usize + 12;
    let w = pow10((k as f64 - e).floor() as usize);
    let m = pow10(k);

    let mut rows = vec![std::iter::once(w.clone()).chain(alphas(k)).collect::<Vec<_>>()];
    for i in 0..n {
        let mut r = vec![Integer::new(); n + 1];
        r[i + 1] = m.clone();
        rows.push(r);
    }
    let mut basis = MpzLatticeBasis::new(rows);

    let _ = mpz_lll_reduce_with_options(&mut basis, MpzLllOptions::default());
    if block >= 2 {
        mpz_bkz(&mut basis, block, 8);
    }

    // Every reduced row with a nonzero first coordinate is a candidate (q, p_1, ..., p_n).
    let mf = m.to_f64();
    let mut best: Option<Approx> = None;
    for row in &basis.basis {
        let q = Integer::from(&row[0] / &w).abs();
        if q == 0 {
            continue;
        }
        let worst = row[1..].iter().map(|x| (x.to_f64() / mf).abs()).fold(0.0, f64::max);
        let cand = Approx { q, worst };
        if best.as_ref().map_or(true, |b| quality(&cand, n) < quality(b, n)) {
            best = Some(cand);
        }
    }
    best.unwrap()
}

const PI: &str = "3.14159265358979323846264338327950288419716939937510582097494459";
const E: &str = "2.71828182845904523536028747135266249775724709369995957496696762";
fn first_primes(n: usize) -> Vec<u32> {
    (2u32..).filter(|&p| (2..p).take_while(|d| d * d <= p).all(|d| p % d != 0)).take(n).collect()
}

fn main() {
    // 1. One number: pi, with q up to about 10^3.
    let a = approximate(1, |k| vec![from_decimal(PI, k)], 3, 0);
    println!("pi: q = {}, |q*pi - p| = {:.3e}", a.q, a.worst);

    // 2. Two and three numbers with one shared denominator, q up to about 10^6.
    let a = approximate(2, |k| vec![from_decimal(PI, k), from_decimal(E, k)], 6, 0);
    println!("(pi, e): q = {}, worst = {:.3e}, quality = {:.3}", a.q, a.worst, quality(&a, 2));
    let a = approximate(3, |k| vec![from_decimal(PI, k), from_decimal(E, k), sqrt_scaled(2, k)], 6, 0);
    println!("(pi, e, sqrt2): q = {}, worst = {:.3e}, quality = {:.3}", a.q, a.worst, quality(&a, 3));

    // 3. Scale up: square roots of the first n primes, Q = 10^n. LLL versus BKZ-20.
    println!("  n  dim   LLL quality   BKZ-20 quality");
    for n in [5usize, 10, 20, 30, 40, 50, 60] {
        let primes = first_primes(n);
        let alphas = |k: usize| primes.iter().map(|&p| sqrt_scaled(p, k)).collect();
        let l = approximate(n, alphas, n as u32, 0);
        let b = approximate(n, alphas, n as u32, 20);
        println!("{n:3}  {:3}   {:11.3}   {:14.3}", n + 1, quality(&l, n), quality(&b, n));
    }
}

Output on an Apple-silicon laptop, in release mode (about three minutes, almost all of it BKZ in the two largest dimensions):

$ cargo run --release
pi: q = 113, |q*pi - p| = 3.014e-5
(pi, e): q = 432134, worst = 3.407e-4, quality = 0.224
(pi, e, sqrt2): q = 760668, worst = 2.063e-3, quality = 0.188
  n  dim   LLL quality   BKZ-20 quality
  5    6         0.346            0.346
 10   11         0.384            0.384
 20   21         0.460            0.499
 30   31         0.545            0.504
 40   41         0.577            0.512
 50   51         0.682            0.633
 60   61         0.710            0.677

One number: recovering 355/113

The first run is the two-dimensional case: $n = 1$ , $Q = 10^3$ , $M = 10^{18}$ , $W = 10^{12}$ . The reduced basis contains the vector

113\,b_0 - 355\,b_1 = (113\cdot 10^{12},\; -30{,}144{,}353{,}364{,}106) \approx (1.13\cdot 10^{14},\; -3.01\cdot 10^{13}),

which says $113\pi - 355 \approx -3.01\cdot 10^{-5}$ . It is the same answer the continued fraction $\pi = [3; 7, 15, 1, 292, \dots]$ gives: the convergents are $3, \tfrac{22}{7}, \tfrac{333}{106}, \tfrac{355}{113}$ , and the next one, $\tfrac{103993}{33102}$ , is beyond $Q = 1000$ . In two dimensions LLL is essentially the Gauss–Lagrange reduction, and reducing this basis retraces the continued fraction algorithm. The quality, $113 \cdot 3.01\cdot 10^{-5} \approx 0.0034$ , is far below 1 because the next partial quotient, 292, is unusually large; $\pi$ is simply very close to this fraction.

Two and three numbers

With $Q = 10^6$ the lattice returns $q = 432134$ for $(\pi, e)$ . Rounding gives $p = (1357589, 1174662)$ , so $\pi \approx 1357589/432134$ and $e \approx 1174662/432134$ with the worse of the two errors at $3.4\cdot 10^{-4}$ (quality 0.224). Adding $\sqrt{2}$ gives $q = 760668$ and $p = (2389709, 2067710, 1075747)$ , worst error $2.1\cdot 10^{-3}$ (quality 0.188). Both qualities sit next to the predicted $0.24$ . We checked these fractions independently with 80-digit decimal arithmetic.

Scaling up to sixty numbers

The last loop approximates $\sqrt{2}, \sqrt{3}, \sqrt{5}, \dots$ , the square roots of the first $n$ primes, with $Q = 10^n$ , and compares LLL alone with LLL followed by BKZ with block size 20:

$n$	dimension	LLL quality	BKZ-20 quality
5	6	0.346	0.346
10	11	0.384	0.384
20	21	0.460	0.499
30	31	0.545	0.504
40	41	0.577	0.512
50	51	0.682	0.633
60	61	0.710	0.677

Three things stand out. First, every entry is below 1: in these dimensions reduction delivers Dirichlet-quality denominators, which a brute-force search could never find at $Q = 10^{60}$ . Second, quality drifts upward with $n$ , as the $\delta^{\,n+1}$ factor predicts. Third, BKZ helps only modestly here and not monotonically (at $n = 20$ the LLL vector happened to be slightly better, since “shorter overall” does not have to mean “smaller maximum error”). BKZ earns its cost in higher dimensions, where LLL’s exponential factor grows. Because this problem has no planted answer, a modest improvement in the root-Hermite factor turns directly into a modest improvement in quality, with no sharp success-or-failure threshold.

Where this is used

Common denominators and rational reconstruction. Measured ratios that should share a small denominator, such as isotope abundances, crystal lattice parameters or recovered floating-point fractions, can be snapped to exact rationals with one shared $q$ .
Gear trains and tunings. Choosing tooth counts that realize several target ratios at once, or an equal temperament whose steps approximate several just intervals (the familiar $\log_2(3/2) \approx 7/12$ generalizes to several intervals at once), is simultaneous approximation.
Resonances and quasi-periodicity. Finding a time $q$ at which several periodic motions nearly realign is the same problem with the frequencies as the $\alpha_i$ .
Theory. Lagarias used exactly this lattice to study the complexity of simultaneous approximation, and it is a standard example of LLL turning a pigeonhole existence argument into an algorithm.

Common mistakes

Forgetting the weight. With $W = 1$ the answer is a valid approximation, but its denominator is whatever the precision $M$ dictates, usually far larger than you want.
Too few digits. If $M$ is not much larger than $Q^{(n+1)/n}$ , truncation noise dominates the error coordinates and the “best” $q$ fits the rounding, not the numbers. Keep guard digits and check the result against higher-precision values.
Using f64 entries. A double holds 15–16 significant digits. Beyond a handful of numbers the basis entries exceed that, and floating-point reduction silently works on a different lattice. Use exact integers.
Reading only the first row. Several reduced rows can have a nonzero first coordinate; scan them all and keep the one with the best maximum error, as the code does.
Expecting uniqueness. Unlike integer relations, there is no single correct $q$ . Different reductions return different, equally valid answers.

Next steps

The same construction, with the unknown integers on one side and scaled residuals on the other, reappears in the integer relation lattice, which transposes the roles of the rows and columns. Adding a target point to be approximated turns it into a closest-vector problem, which is the subject of CVP by Kannan embedding. For why a short vector must exist at all, see Minkowski’s theorem in the geometry of numbers.

7.3 Simultaneous Diophantine approximation