A minimal implementation of the F4 algorithm for Gröbner bases, over the
prime field
Given -F QQ the same computation runs modulo a
sequence of primes and the answer is reconstructed over
Polynomial arithmetic comes from FLINT
(nmod_mpoly / fmpq_mpoly); the sparse elimination comes from
SparseRREF. Everything else — selection,
symbolic preprocessing, the Macaulay matrix, the structured elimination, the
criteria, the extraction, the multi-modular driver and the FLINT oracle — is in
f4.hpp.
--verify cross-checks the result against FLINT — the computed basis is
handed to FLINT's own fmpz_mod_mpoly machinery, which re-checks the
S-polynomials, the containment of the input and is_groebner. It is off by
default, because on a large input the check costs more than the computation.
-
One header. Drop
f4.hppinto a project; the only dependencies are FLINT and the header-only SparseRREF. A ~570-linemain.cppis a CLI around it, not a requirement. -
Independently verifiable.
--verifyhands the basis to FLINT's ownfmpz_mod_mpolymachinery, which re-checks the S-polynomials, the ideal containment andis_groebnerin a separate context. - Three elimination kernels picked per batch from the matrix shape.
-
$\mathbb{F}_p$ and$\mathbb{Q}$ through the same code path — over$\mathbb{Q}$ the whole F4 is re-run per prime and only the final basis is reconstructed.
All timings below were measured on one Windows machine (MinGW-w64 g++ 16.1.0,
vcpkg FLINT 3.6.0, -t 4, p = 32003). "Ours" is the time the F4 loop reports
and excludes process startup; the oracle runs only with --verify. Singular
4.4.1 runs under WSL and its time is its own std call (measured inside
Singular, whose rtimer is raised to 1 µs resolution); Mathematica 15 is
GroebnerBasis with AbsoluteTiming.
| system | vars | GB | ours | Singular | Mathematica |
|---|---|---|---|---|---|
| katsura-3 | 4 | 7 | 0.002 s | 0.000 s | 0.001 s |
| katsura-4 | 5 | 13 | 0.003 s | 0.000 s | 0.004 s |
| katsura-5 | 6 | 21 | 0.005 s | 0.001 s | 0.018 s |
| katsura-6 | 7 | 38 | 0.014 s | 0.007 s | 0.189 s |
| cyclic-4 | 4 | 7 | 0.003 s | 0.000 s | 0.003 s |
| cyclic-5 | 5 | 20 | 0.012 s | 0.001 s | 0.009 s |
| cyclic-6 | 6 | 45 | 0.026 s | 0.005 s | 0.151 s |
| cyclic-7 | 7 | 209 | 0.909 s | 0.631 s | 20.38 s |
| noon-5 | 5 | 72 | 0.103 s | 0.003 s | 0.030 s |
| noon-6 | 6 | 187 | 2.331 s | 0.023 s | 0.420 s |
| reimer-4 | 4 | 17 | 0.009 s | 0.001 s | 0.013 s |
| reimer-5 | 5 | 38 | 0.065 s | 0.026 s | 0.286 s |
| reimer-6 | 6 | 95 | 1.472 s | 5.076 s | 18.60 s |
13/13 agree with Singular (each basis reduces to 0 modulo the other, decided by Singular) and 13/13 match Mathematica element by element.
Two things about these matrices explain most of the design. cyclic-7's are
tall (2.2× as many rows as columns) and mostly leading columns — only 433 of
those 2 985 columns are non-pivot columns (at most 498 over all 33 batches). The
default kernel keeps a dense accumulator exactly that wide, which is why it wins
there. The noon systems are the opposite shape (noon-6: 50 122 × 50 629, half
a million nonzeros) and they go to SparseRREF — and they are also where
Singular beats this engine by two orders of magnitude: its Buchberger runs
in 23 ms where a batch F4 needs 2.3 s. That is the honest shape of the
comparison: F4 wins big on the tall, over-constrained, mostly-leading-column
matrices it was built for, and loses on square dense ones.
| system | GB | primes | modulus | ours | Singular | Mathematica |
|---|---|---|---|---|---|---|
| katsura-3 | 7 | 1 | 61 bits | 0.002 s | 0.000 s | 0.00 s |
| katsura-4 | 13 | 2 | 121 bits | 0.006 s | 0.001 s | 0.00 s |
| katsura-5 | 21 | 5 | 301 bits | 0.022 s | 0.005 s | 0.01 s |
| katsura-6 | 38 | 14 | 841 bits | 0.160 s | 0.104 s | 0.16 s |
| cyclic-4 | 7 | 1 | 61 bits | 0.003 s | 0.000 s | 0.00 s |
| cyclic-5 | 20 | 1 | 61 bits | 0.011 s | 0.001 s | 0.01 s |
| cyclic-6 | 45 | 2 | 121 bits | 0.050 s | 0.029 s | 0.17 s |
| cyclic-7 | 209 | 10 | 601 bits | 10.7 s | > 4 min | not measured |
| noon-5 | 72 | 1 | 61 bits | 0.132 s | 0.003 s | 0.03 s |
| reimer-5 | 38 | 2 | 121 bits | 0.137 s | 0.108 s | 0.40 s |
9/9 agree with Singular. cyclic-7 is the ceiling here: Singular did not
finish it in four minutes and the Mathematica reference was never computed —
over
The number of primes is decided by the height of the final basis only, and
reconstruction itself is negligible (0.14 s for all 27 187 coefficients of
cyclic-7). Over
The primes themselves are as wide as the engine allows: --qbits (default 60)
sets the width of the first one, and the only constraint is SparseRREF's
Requirements:
- a C++20 compiler;
- FLINT 3.x with GMP and MPFR;
- SparseRREF, header-only, at
./SparseRREF.
f4.exe [options] [generator ...]
Generators are parsed by FLINT's nmod_mpoly_set_str_pretty (or
fmpq_mpoly_set_str_pretty with -F QQ), e.g. "x^2+y^2+z^2-1". With no
generator and no -e, the help is printed.
Basics
| option | meaning |
|---|---|
-p, --prime P |
prime modulus, below 2⁶³ (default 32003) |
-v, --vars LIST |
comma-separated variable names (default: guessed from the input) |
-g, --gen POLY |
a generator; may be repeated (same as a positional argument) |
--file F |
read generators from a file, one polynomial per line (# comments; F = - reads stdin). A bare argument naming a file is read the same way |
-e, --example N |
run built-in example N (0..3) |
-t, --threads N |
threads used by the elimination and the matrix fill-in (default 1) |
--qbits N |
width in bits of the first prime the |
-F, --field F |
coefficient field: Zp (default) or QQ
|
-h, --help |
the option list |
Verification and output
| option | meaning |
|---|---|
--verify |
cross-check the result against FLINT (off by default) |
-V, --verbose |
print FLINT's reduced GB when the two disagree (implies --oracle) |
--oracle |
additionally check that FLINT's own naive Buchberger gives the same reduced basis, element by element; the slow part on large inputs |
--dump-gb |
print the basis as machine-readable POLY<TAB>... lines |
Elimination kernels (the default is shape-aware)
| option | meaning |
|---|---|
--sparserref |
force the general sparse RREF |
--gbla-style |
force the GBLA-style structured elimination (the default kernel) |
--buckets |
structural bucket sweep: forward scan only, fast on small matrices, much slower on large ones |
--backsub |
also do the back substitution (full RREF); only SparseRREF supports it, and selecting it switches there |
--all-divisors |
one closure row per divisor (textbook F4; ~24× bigger matrices, kept for comparison) |
Criteria
| option | meaning |
|---|---|
--lazycrit / --nolazycrit |
re-apply the chain criterion when a batch is selected (on by default) |
--bguard / --nobguard |
classical strictness conditions on the backward criterion (off by default) |
Examples
.\f4.exe # print the help
.\f4.exe -e 1 # cyclic-3
.\f4.exe -e 1 --verify # ... checked against FLINT
.\f4.exe "x^2+y^2+z^2-1" "x*y-1" # guess the variables
.\f4.exe -v x,y,z,w -p 65521 "x+y+z+w" "x*y+y*z+z*w+w*x" "x*y*z*w-1"
.\f4.exe -F QQ "x^2+y^2-1" "x*y-1" # multi-modular over Q
.\f4.exe -F QQ "x^3+y^3-1" "x*y-2" # a rational answerEach batch also prints a statistics line (RREF, LA, GBLA, BATCH) with
the matrix shape, the kernel that was chosen and where the time went. There is
no matrix export by design.
Using the header directly:
#include "f4.hpp"
nmod_mpoly_ctx_t ctx;
nmod_mpoly_ctx_init(ctx, 3, ORD_DEGREVLEX, 32003);
f4::F4 engine(ctx, 32003);
engine.set_threads(4);
std::vector<f4::Poly> G = engine.run(generators); // reduced GBOne matrix elimination per degree level, exactly as in Faugère's F4:
G <- input, made monic (no initial minimisation -- see below)
pairs <- critical pairs of G, filtered by the Gebauer-Moeller criteria
while pairs is not empty:
d <- minimal lcm degree among pairs
P <- the pairs with lcm degree d
P <- P filtered *again* by the chain criterion (lazy criteria)
# ---- symbolic preprocessing ----
rows <- { (lcm/L_i) g_i , (lcm/L_j) g_j : (i,j) in P }
M <- every monomial occurring in those rows, bucketed by total degree
sweep M from high degree to low:
for m in M, for g in G with L_g | m:
rows <- rows + { (m/L_g) g }; break # only the *first* divisor
cols <- M sorted by degrevlex, descending
# ---- linear algebra ----
A <- Macaulay matrix (rows x cols) over F_p
kernel<- SparseRREF if npc > 1024 or npc*(rows-npiv) > 3*nnz
GBLA-style otherwise
A <- echelon(A)
if a nonzero constant appears: return {1} # the ideal is the whole ring
# ---- extraction ----
G <- G + { nonzero rows of A whose leading monomial
is not divisible by any L_g of G }
pairs <- pairs + the new elements' pairs, GM-filtered
return interreduce(G)
The final interreduction (minimise, then reduce the tails) returns the unique
reduced Gröbner basis. Rows are extracted from an echelon form rather than a
full RREF on purpose: the extraction rule reads only "did this row become zero"
and "what is its leading monomial", and both are already final in the echelon
form. --backsub is kept for comparison.
Three kernels, all exact over f4.hpp:
-
GBLA-style structured elimination (default). The matrix is cut by pivot
columns into A/B/C/D blocks: A is turned into the identity, one pass of
[I | B']zeroes C, and the remaining dense block D goes to FLINT's densenmod_mat_rref. The accumulator that zeroes C is as wide as the number of non-pivot columns and uses delayed modular reduction:ulongwhen$(ncol+1)(p-1)^2 < 2^{64}$ , otherwiseunsigned __int128with a reduction every$\lfloor max/(p-1)^2 \rfloor$ terms; where the compiler has no 128-bit integer at all (MSVC), that last case reduces after every term instead, which gives the same answer with more reductions per row. The algorithm is the one described by Faugère–Lachartre's GBLA; no GBLA code is used and no data is exchanged with it. - SparseRREF — the general fill-in-aware sparse RREF, used when the shape test says the dense accumulator would lose.
-
Bucket sweep (
--buckets) — a forward-scan-only reference implementation. It wins on small matrices (15× on cyclic-4) and loses badly on large ones (5× slower on cyclic-6, unfinished on cyclic-7): on a 717-row matrix the cascade still performs 206 927 row eliminations.
The dispatch is per batch (npc > 1024 or npc·(rows − npiv) > 3·nnz →
SparseRREF). npc·(rows − npiv) is the number of dense cells the accumulator
kernel sweeps; SparseRREF's work follows the actual nonzeros, so the ratio of
the two is the shape test. Both paths produce an echelon form, so mixing them
across batches is sound.
Symbolic preprocessing holds most of the remaining constant factor:
g_terms, a cache of each basis element's terms as (exponent delta from the leading term, coefficient), so a row's monomials are vector additions instead of FLINT unpacking. It carries coefficients, so it is rebuilt per prime.MonoSet/MonoMap, open-addressed monomial sets whose table stores{32-bit hash, 32-bit index}into a flat vector of keys: one probe touches 8 contiguous bytes instead of a 74-byte key, and full comparison happens only on a real hash collision.- the closure swept by total-degree buckets instead of a degrevlex heap;
LFilt, a packed{mask, degree}filter array for the inner divisor scan.
The engine's arithmetic is fmpq_reconstruct_fmpz, accepting only reconstructions with leading
coefficient 1). A prime whose basis has a different support is discarded as
unlucky, and a fresh unused prime re-checks the reconstructed basis. Entry
point: F4::run_over_q; the check is F4::verify_gb_over_q.
The oracle that ships here (--verify). verify_mod_p (at the end of f4.hpp)
converts the computed basis into FLINT's fmpz_mod_mpoly in a fresh context —
sharing nothing with the engine except the string spelling — and checks three
things: every S-pair of the basis reduces to zero, every input generator reduces
to zero, and FLINT's own is_groebner agrees. Points 1 and 2 together are a
proof that the output is a Gröbner basis of the input ideal. --oracle
additionally compares against FLINT's naive Buchberger element by element.
External cross-checks (not part of this repository): every
system in the tables above was also run past Singular's std over
GroebnerBasis, whose result must agree element by element.
Neither tool shares any code with this engine.
Three lessons from that harness that shaped the defaults:
- an optimisation that deleted redundant basis elements inside the main loop was 2× faster on cyclic-6 with every benchmark green, but failed 1–3 % of random systems — it is unsound, and it was removed;
- a benchmark driver bug sent options to only the first instance, which produced two false negative results before it was noticed;
- the delayed-reduction accumulator used to flush on
++terms == limit, while its term counter is seeded with 1 when a sweep starts from a carried-in partial. For$p > 2^{63.5}$ the limit is 1, so the test never fired and a threaded run silently overflowed — returning a spurious{1}that the FLINT oracle cannot refute, because with$G = {1}$ the containment check is vacuous. Adding the portable (no-__int128) path for MSVC is what exposed it: the two paths disagreed, and single-threaded runs plus every prime below$2^{63.5}$ were fine. The flush test is now>=, and the two paths are differentially tested against each other.
| file | contents |
|---|---|
f4.hpp |
the whole library: monomials, symbolic preprocessing (MonoSet, LFilt, degree buckets), the three elimination kernels, the F4 engine, the Gebauer–Möller criteria, the multi-modular verify_mod_p
|
main.cpp |
CLI and the built-in examples |
Nothing needs configuring. The modulus has to be a prime below n_mulmod_shoup multiplication needs it) and the
engine inherits it. The only dependency beyond FLINT is SparseRREF; the engine's
own kernels do not need it, but the default dispatch calls into it.
-
Faugère's
$F$ criterion / signature-based algorithms. Deliberately omitted rather than risk a silently wrong criterion. -
A rational-function coefficient field. The algorithm is field-agnostic;
the practical route would be evaluation and interpolation over
$\mathbb{F}_p$ rather than a new scalar type.
MIT — see LICENSE.
- J.-C. Faugère, A new efficient algorithm for computing Gröbner bases (F4), J. Pure Appl. Algebra 139 (1999) 61–88.
- SparseRREF — header-only sparse exact RREF, MIT. The linear algebra backbone.
- FLINT — Fast Library for Number Theory, LGPL.
nmod_mpoly/fmpq_mpolyarithmetic and the independent oracle. - GBLA — Faugère–Lachartre's structured elimination (LGPL 2.1). Read for the algorithm only; no code or data is shared.