7.2 Integer relations: finding a minimal polynomial with LLL

Suppose you computed $\alpha \approx 3.14626436994197$ and suspect it is algebraic. Is there a polynomial with small integer coefficients that vanishes at $\alpha$? This is an integer relation problem, and it is the first embedding in this series where a single choice, the scaling factor $N$, decides between the right answer and garbage. We follow the recipe from 7.1 Embedding basics and solve it with LLL from g6k-rs.

The problem

Given real numbers $x_0, \dots, x_d$, find integers $c_0, \dots, c_d$, not all zero and preferably small, with

For a minimal polynomial we take $x_i = \alpha^i$. A relation is then a polynomial $p(X) = \sum_i c_i X^i$ with $p(\alpha) = 0$. We only know the $x_i$ to finite precision, so the relation holds approximately in our data and exactly for the true numbers.

Building the lattice

Step 1 of the recipe asks for the small object: the coefficient vector $c$. Step 2 writes the relation as a linear form, $\sum_i c_i x_i$, which should be zero. Real numbers cannot go into an integer lattice, so we scale them by a large integer $N$ and round. Using the row convention (rows of $B$ are basis vectors), row $i$ is the unit vector $e_i$ followed by one extra coordinate:

The integer combination with coefficients $c$ is the lattice vector

For the true relation, $\sum_i c_i x_i = 0$, so $r$ is only the accumulated rounding error: $|r| \le \tfrac12 \sum_i |c_i|$ when the $x_i$ are exact. A vector $c$ that is not a relation has $r \approx N \sum_i c_i x_i$, which is large. The last column acts as a penalty that makes every non-relation long. The target norm is

How large must N be?

The lattice has rank $n = d+1$ in ambient dimension $d+2$. It is not full rank, so its determinant is $\sqrt{\det(BB^\top)}$. The Gram matrix is the identity plus a rank-one term, which gives

Step 5 of the recipe compares the target with the Gaussian heuristic:

Read the right side as a rule of thumb: the number of digits of N grows linearly with the degree, roughly $n \log_{10}(t\sqrt{2\pi e/n})$. Doubling the degree roughly doubles the digits you need. This is a heuristic for random-looking lattices, not a proof, and the examples below show it is a useful but not sharp guide.

Example 1: √2 + √3, degree 4

The number $\alpha = \sqrt2 + \sqrt3$ satisfies $\alpha^4 - 10\alpha^2 + 1 = 0$, so the target is $c = (1, 0, -10, 0, 1)$ with $\|c\| = \sqrt{102} \approx 10.1$. Here $n = 5$ and $\|x\| = \|(1, \alpha, \dots, \alpha^4)\| \approx 103.3$, so the rule above gives $N \gtrsim 2.2 \times 10^4$. The program below builds the basis in f64 and runs LLL for several values of $N$. Add this to Cargo.toml:

[dependencies]
g6k-rs = "=0.1.0"
# only needed for the high-precision example
rug = { version = "1.26", default-features = false, features = ["integer", "float"] }

// Recover the minimal polynomial of alpha = sqrt(2) + sqrt(3) with LLL.
// Rows of the basis are lattice vectors: row i = (e_i | round(N * alpha^i)).
use g6k_rs::{LLLParams, LatticeBasis, lll_reduce};

fn main() {
    let alpha = 2f64.sqrt() + 3f64.sqrt();
    let degree = 4;
    let dim = degree + 1;
    for exp in [2, 4, 6, 8] {
        let n_scale = 10f64.powi(exp);
        let rows: Vec<Vec<i64>> = (0..dim)
            .map(|i| {
                let mut row = vec![0i64; dim + 1];
                row[i] = 1;
                row[dim] = (n_scale * alpha.powi(i as i32)).round() as i64;
                row
            })
            .collect();
        let mut basis = LatticeBasis::from_integers(&rows);
        let result = lll_reduce(&mut basis, &LLLParams::default());
        // The first reduced row is the shortest vector found: (c_0, ..., c_4, residual).
        let v: Vec<i64> = basis.basis[0].iter().map(|&x| x.round() as i64).collect();
        println!(
            "N = 1e{exp}  converged = {:<5}  c = {:?}  residual = {}",
            result.converged, &v[..dim], v[dim]
        );
    }
}

Output with g6k-rs 0.1.0:

N = 1e2  converged = true   c = [-1, 2, -1, -3, 1]  residual = -3
N = 1e4  converged = true   c = [1, 0, -10, 0, 1]  residual = -2
N = 1e6  converged = true   c = [1, 0, -10, 0, 1]  residual = 5
N = 1e8  converged = false  c = [1, 0, 0, 0, 0]  residual = 100000000

The table adds $N = 10^{3.5}$ from the same program run on a finer grid (10f64.powf(exp)), with the heuristic values computed from the formula above.

Three lessons are in that table. For small $N$ the penalty column is too cheap, and LLL correctly returns a vector shorter than our target: $-1 + 2\alpha - \alpha^2 - 3\alpha^3 + \alpha^4 \approx -0.051$ is small but not zero, and at $N = 100$ that is good enough to look like a relation. Near the heuristic threshold the answer appears (here already at $10^4$, a bit below the predicted $2.2 \times 10^4$, which is normal for a heuristic). And at $N = 10^8$ the floating-point reducer gives up.

Why f64 stops working

The floating-point LatticeBasis stores entries as f64 and computes the Gram–Schmidt data in double precision. With $N = 10^8$ the last column holds values near $10^{10}$, squared lengths near $10^{20}$, while the Gram–Schmidt lengths we care about are about 1. Double precision keeps only about 16 significant digits, so those small quantities vanish in cancellation. g6k-rs detects the non-positive Gram–Schmidt length and returns the basis unreduced with converged = false instead of a wrong answer. Always check that flag.

With f64 the working window for this example is roughly $10^4 \le N \le 10^6$. That window closes as the degree grows, because the lower end rises with the degree and the upper end is fixed by double precision. For anything beyond a toy, use exact integers.

Example 2: 2^(1/3) + 3^(1/3), degree 9, with exact integers

Let $\alpha = \sqrt[3]{2} + \sqrt[3]{3} \approx 2.70217$. Cubing gives $\alpha^3 = 5 + 3\sqrt[3]{6}\,\alpha$, and eliminating the cube root gives $(\alpha^3 - 5)^3 = 162\,\alpha^3$, that is

The target $c = (-125, 0, 0, -87, 0, 0, -15, 0, 0, 1)$ has norm $\sqrt{23420} \approx 153$, the rank is $n = 10$, and $\|x\| \approx 8268$. The heuristic now asks for $N \gtrsim 1.2 \times 10^{19}$, and the largest entry $N\alpha^9$ is around $10^{23}$. That no longer fits in an i64 (maximum about $9.2 \times 10^{18}$) and is far beyond the $2^{53} \approx 9 \times 10^{15}$ integers that f64 represents exactly. We switch to MpzLatticeBasis, which stores GMP integers and runs LLL with exact row operations, and compute $\alpha$ with MPFR at 256 bits.

// Degree-9 integer relation with exact integers: alpha = 2^(1/3) + 3^(1/3).
// The scaled column needs about 24 digits, beyond i64 and f64, so use MpzLatticeBasis.
use g6k_rs::{MpzLatticeBasis, MpzLllOptions, mpz_lll_reduce_with_options};
use rug::ops::Pow;
use rug::{Float, Integer};

fn main() {
    let prec = 256; // bits of working precision for alpha, about 77 decimal digits
    let alpha = Float::with_val(prec, 2u32).cbrt() + Float::with_val(prec, 3u32).cbrt();

    let degree = 9;
    let dim = degree + 1;
    for digits in [15u32, 18, 20, 22, 24] {
        let n_scale = Float::with_val(prec, Integer::from(10u32).pow(digits));
        let basis_rows: Vec<Vec<Integer>> = (0..dim)
            .map(|i| {
                let mut row = vec![Integer::new(); dim + 1];
                row[i] = Integer::from(1);
                let scaled = Float::with_val(prec, &n_scale * alpha.clone().pow(i as u32));
                row[dim] = scaled.round().to_integer().unwrap();
                row
            })
            .collect();
        let mut basis = MpzLatticeBasis { rows: dim, cols: dim + 1, basis: basis_rows };
        let result = mpz_lll_reduce_with_options(&mut basis, MpzLllOptions::default());
        let v = &basis.basis[0];
        println!(
            "N = 1e{digits}  converged = {}  c = {:?}  residual = {}",
            result.converged,
            v[..dim].iter().map(|x| x.to_i64().unwrap_or(i64::MAX)).collect::<Vec<_>>(),
            v[dim]
        );
    }
}

Output with g6k-rs 0.1.0:

N = 1e15  converged = true  c = [-27, -19, 37, -11, -25, 1, 4, -45, -24, 15]  residual = -18
N = 1e18  converged = true  c = [-125, 0, 0, -87, 0, 0, -15, 0, 0, 1]  residual = 23
N = 1e20  converged = true  c = [-125, 0, 0, -87, 0, 0, -15, 0, 0, 1]  residual = 23
N = 1e22  converged = true  c = [-125, 0, 0, -87, 0, 0, -15, 0, 0, 1]  residual = 25
N = 1e24  converged = true  c = [-125, 0, 0, -87, 0, 0, -15, 0, 0, 1]  residual = 4

At $N = 10^{15}$ (planting ratio $\mathrm{gh}/t \approx 0.39$) LLL returns a shorter spurious relation. From $10^{18}$ on (ratio about 0.78 and rising) the minimal polynomial appears. Again the true threshold sits a little below the heuristic one, and again the safe choice is a margin above it.

When you only know some digits

In the program above $\alpha$ is known to 77 digits. In practice the input is a measured or computed value with limited accuracy. Truncating $\alpha$ to $k$ decimal digits and fixing $N = 10^{22}$ gives:

This matches the digits rule: $\log_{10} N = 22$ plus $\log_{10}(9\,\alpha^8) \approx 4.4$ suggests about 26 digits, and the relation appears from 24. With fewer digits, rounding noise in the last column is large enough that LLL finds vectors that fit the noise.

Checking the answer

LLL always returns some short vector, so a result is a candidate, not a proof. Evaluate it at much higher precision than you used to find it. A genuine relation drops to the new precision limit; a spurious one stays at the level of the old noise:

from decimal import Decimal as D, getcontext
getcontext().prec = 60
alpha = D(2) ** (D(1) / 3) + D(3) ** (D(1) / 3)

for c in ([-125, 0, 0, -87, 0, 0, -15, 0, 0, 1],       # found with 24 known digits
          [1, -145, 112, -7, -196, 204, -37, 39, 146, -60]):  # found with 22 known digits
    print(c, "%.3E" % sum(D(ci) * alpha**i for i, ci in enumerate(c)))

# [-125, 0, 0, -87, 0, 0, -15, 0, 0, 1] -2.000E-56
# [1, -145, 112, -7, -196, 204, -37, 39, 146, -60] -9.411E-18

Two further checks turn a candidate into a statement about the minimal polynomial. The candidate must be irreducible over the integers (any computer algebra system can factor it), and you need an independent argument that $\alpha$ is algebraic of at most that degree, like the cube-root elimination above. The lattice alone only says that no relation shorter than a certain bound was found.

PSLQ, the classical alternative

The PSLQ algorithm of Ferguson and Bailey solves the same problem without choosing $N$: it works directly with the real numbers and, run with enough precision, finds a relation or proves a lower bound on the norm of any relation. It is the standard tool in experimental mathematics (mpmath ships pslq and findpoly). The lattice approach is worth knowing because it generalizes: the same penalty-column construction handles several equations at once, weighted coordinates, and extra congruence constraints, which is where the next guides go.

Common mistakes

• N too small. The penalty column becomes cheap, and LLL returns a genuine shorter vector that is not a relation. Use the heuristic threshold plus a margin.
• N larger than your digits support. Scaling by $10^{22}$ numbers that are only correct to 18 digits multiplies noise, not information.
• Floating-point overflow of precision. f64 bases break once squared entries exceed about $10^{16}$ times the small Gram–Schmidt lengths. Check converged and move to MpzLatticeBasis.
• Forgetting the constant term. Include $x_0 = 1$. Without it, only relations with $c_0 = 0$ can be found, and the minimal polynomial of $\sqrt2 + \sqrt3$ has $c_0 = 1$.
• Trusting the first row blindly. Verify at higher precision and check irreducibility before calling a result the minimal polynomial.

Next, 7.3 Simultaneous Diophantine approximation turns the construction around: instead of one equation among many numbers, it asks for one denominator that fits many numbers at once, and introduces coordinate weights.

Further reading

Lenstra, Lenstra and Lovász, Factoring polynomials with rational coefficients (1982), where LLL is introduced for exactly this kind of problem. Kannan, Lenstra and Lovász, Polynomial factorization and nonrandomness of bits of algebraic and some transcendental numbers (1988), for the precision bounds. Ferguson, Bailey and Arno, Analysis of PSLQ, an integer relation finding algorithm (1999). For the reduction side, see the LLL chapter.