AutoNumerics: Automated Numerical Design

• AutoNumerics is a research motif that automates numerical design by generating, verifying, and executing constructs traditionally handcrafted by experts.
• It encompasses diverse methodologies including data-driven invariant synthesis, precision tracking for floating-point operations, evolutionary search for function approximations, and autonomous PDE solver pipelines.
• These automated approaches yield practical benefits in computational accuracy and efficiency, as shown by rigorous benchmarks and innovative code-generation techniques.

to=arxiv_search üpjil {"query":"AutoNumerics", "max_results": 10} AutoNumerics is a label used in the arXiv literature for several distinct lines of research that automate parts of numerical reasoning, numerical software construction, and numerical representation. The term has been applied to data-driven numerical invariant synthesis for program verification, precision tracking for IEEE-754 computations, automated discovery of transcendental-function approximations, and autonomous generation of PDE solvers from natural-language specifications (Bouajjani et al., 2022). Related antecedents include automatic synthesis of abstract transformers for template numerical constraints and semiautomatic ring constructions used as “automatic numeric systems” for geometry (0909.4013). Taken together, these works do not define a single canonical framework; rather, they mark a family of efforts to replace hand-crafted numerical design with algorithmic generation, verification, and execution.

1. Terminological scope and research landscape

In the recent literature, the name AutoNumerics has been attached to at least four substantially different technical agendas. One line studies numerical invariant synthesis from ICE samples using decision trees and automatically generated attributes (Bouajjani et al., 2022). A second line introduces algorithms and data structures for numerical computations with automatic precision estimation, centered on carrying an auxiliary correctness counter alongside every floating-point value (Netay, 2024). A third line, "AutoNumerics-Zero," uses evolutionary search over compute graphs to discover approximation programs for transcendental functions such as $2^x$ under fixed precision or speed objectives (Real et al., 2023). A fourth line presents an autonomous, PDE-agnostic multi-agent pipeline that designs, implements, debugs, and verifies classical numerical solvers for PDEs directly from natural-language descriptions (Du et al., 19 Feb 2026).

Older adjacent work broadens the conceptual background. Monniaux’s method for "automatic modular abstractions for template numerical constraints" automatically generates optimal abstract transformers and least inductive invariants in template linear-constraint domains by quantifier elimination (0909.4013). A different strand studies semiautomatic rings as dense subrings of $\mathbb{R}$ whose addition and comparisons are automatic and whose multiplication is semiautomatic, explicitly describing them as an “automatic numeric system” for geometric constructions (Gao et al., 2021).

A plausible implication is that “AutoNumerics” should be understood as a research motif rather than a uniquely identified method. Across these works, the unifying theme is automation of numerical artifacts that have traditionally required manual mathematical design: attributes for invariants, precision bookkeeping, transcendental kernels, abstract transformers, or PDE discretizations.

Usage	Core object	Primary mechanism
Numerical invariant synthesis	Program invariants	ICE-DT with separator-derived attributes
Automatic precision estimation	Floating-point computations	Per-value exact_bits counters
AutoNumerics-Zero	Mathematical function approximations	Evolutionary search over compute graphs
PDE solver generation	Classical PDE solvers	Multi-agent LLM pipeline
Automatic modular abstractions	Abstract transformers and invariants	Quantifier elimination in LRA
Automatic numeric systems	Dense semiautomatic rings	Finite-automaton arithmetic

2. Data-driven numerical invariant synthesis

In "Data-driven Numerical Invariant Synthesis with Automatic Generation of Attributes," AutoNumerics is a learning-and-verification procedure for numerical invariants based on the ICE-DT schema (Bouajjani et al., 2022). The input is an ICE-sample

$S=(S^+,S^-,S^\to),$

where $S^+\subseteq\mathbb{Z}^n$ contains positive examples, $S^-\subseteq\mathbb{Z}^n$ contains negative examples, and $S^\to$ contains implications $p\to q$ induced by program transitions. Consistency requires

$$\forall(p\to q)\in S^\to:\ (p\in S^+ \Rightarrow q\in S^+) \ \text{and}\ (q\in S^- \Rightarrow p\in S^-).$$

A decision tree over a finite attribute set $A$ is considered consistent with the sample when it labels all positive states with $+$, all negative states with $\mathbb{R}$0, and respects the implication condition $\mathbb{R}$1.

The distinctive contribution is automatic attribute generation. Rather than relying on blind template enumeration, the method constructs an $\mathbb{R}$2-separator in a chosen numerical abstract domain

$\mathbb{R}$3

where $\mathbb{R}$4 may be Intervals, Octagons, or Polyhedra (Bouajjani et al., 2022). A separator is a finite family $\mathbb{R}$5 that covers positive states, excludes negative ones, and is consistent with implications. The construction begins from the basic separator $\mathbb{R}$6 and repeatedly joins elements $\mathbb{R}$7 whenever $\mathbb{R}$8, adding singleton implication targets as needed. The process stops when no further join is possible, yielding a join-maximal separator.

The atomic linear constraints defining each separator element become candidate attributes:

$\mathbb{R}$9

These attributes are then supplied to ICE-DT, whose decision-tree generalization can discard some constraints on some paths and thereby move from bounded convex separator components to potentially unbounded invariants (Bouajjani et al., 2022). The outer synthesis loop alternates between learning a candidate invariant $S=(S^+,S^-,S^\to),$0, checking inductiveness via SMT on

$S=(S^+,S^-,S^\to),$1

and updating the sample with positive, negative, or implication counterexamples.

The implementation, "NIS" (Numerical Invariant Synthesizer), is written in C++, uses the APRON library for Intervals, Octagons or Polyhedra, and uses Z3 as SMT-solver (Bouajjani et al., 2022). On 164 LIA benchmarks from SyGuS-Comp’19 with timeout $S=(S^+,S^-,S^\to),$2 s, NIS(poly) solves 154/164, and the virtual best of {NIS(int),NIS(oct),NIS(poly)} solves 160/164; the comparison figures reported are LoopInvGen solves 138, CVC4 solves 129, Spacer solves 136, and original ICE-DT solves 122 (Bouajjani et al., 2022). Runtime profiles are reported to show that NIS(poly) is often faster than LoopInvGen and CVC4, and competitive with Spacer.

This formulation places AutoNumerics within program verification rather than numerical analysis in the narrow sense. Its central problem is not solving a PDE or approximating a function, but inferring inductive numerical assertions from semantic samples. The method’s dependence on a chosen convex abstract domain suggests that its automation is data-guided rather than domain-free.

3. Automatic precision estimation for IEEE-754 computations

A different work uses AutoNumerics to denote a numerical-accuracy tracking scheme for floating-point pipelines (Netay, 2024). Its central representation extends each floating-point variable $S=(S^+,S^-,S^\to),$3 to $\mathbb{R}$37 where exact_bits = m means that “the top $S=(S^+,S^-,S^\to),$4 bits of value’s 52-bit significand agree with the ideal real result” (Netay, 2024). The auxiliary counter is denoted $S=(S^+,S^-,S^\to),$5 and records “how many mantissa bits are still guaranteed correct.”

For addition and subtraction, the method models both mantissa alignment and carry propagation. If $S=(S^+,S^-,S^\to),$6, with alignment shift

$S=(S^+,S^-,S^\to),$7

then conservatively

$S=(S^+,S^-,S^\to),$8

after which one subtracts one more bit for the selected rounding mode (Netay, 2024). For multiplication,

$S=(S^+,S^-,S^\to),$9

with the final “$S^+\subseteq\mathbb{Z}^n$0” accounting for rounding. Division is treated analogously.

Elementary functions are handled through condition numbers. If $S^+\subseteq\mathbb{Z}^n$1 is differentiable at $S^+\subseteq\mathbb{Z}^n$2, the paper defines

$S^+\subseteq\mathbb{Z}^n$3

and assigns

$S^+\subseteq\mathbb{Z}^n$4

All transcendental calls are treated as black boxes whose $S^+\subseteq\mathbb{Z}^n$5 is assumed known in closed form (Netay, 2024). Rounding is analyzed by comparing the first flipped bit with the first inexact bit already present. The framework also specifies rules for max(x,y), pairwise tree sums over arrays, and matrix multiplication.

The matrix-multiplication component is notable because it avoids re-execution in “exactness mode.” If $S^+\subseteq\mathbb{Z}^n$6 carry per-entry counters $S^+\subseteq\mathbb{Z}^n$7, then the fast “tropical” bound is

$S^+\subseteq\mathbb{Z}^n$8

with a stronger but costlier alternative

$S^+\subseteq\mathbb{Z}^n$9

Theorem 5.2 is reported to refine this through

$S^-\subseteq\mathbb{Z}^n$0

Backpropagation reuses the same matrix-multiplication bounds, and convolutions are reduced to im2col + matrix multiply (Netay, 2024).

The implementation trade-offs are explicit. The reported storage overhead is “one extra byte per float $S^-\subseteq\mathbb{Z}^n$1,” pointwise arithmetic is “~20–30 % slower,” and reduction/MatMul overhead is “1–3 extra BLAS calls on integer/weighted data vs your original call” (Netay, 2024). In a GPU context, this can be “10–15 % extra wall-time.” The method also proposes a safety offset $S^-\subseteq\mathbb{Z}^n$2 to accommodate nondeterminism, noting that “in practice $S^-\subseteq\mathbb{Z}^n$3 or $S^-\subseteq\mathbb{Z}^n$4 bits is enough to avoid false ‘exact’ reports on OpenMP or CUDA races.”

This version of AutoNumerics is an online error-accounting discipline rather than an overview framework. It does not search for better numerical algorithms; instead, it instruments existing computations so that catastrophic cancellation and precision loss can be detected without exact recomputation.

4. AutoNumerics-Zero and the automated discovery of function approximations

"AutoNumerics-Zero: Automated Discovery of State-of-the-Art Mathematical Functions" recasts the term as a framework for automatically discovering approximation programs for transcendental functions from scratch (Real et al., 2023). The target function described in detail is

$S^-\subseteq\mathbb{Z}^n$5

on either the real interval $S^-\subseteq\mathbb{Z}^n$6 or the float32 values in $S^-\subseteq\mathbb{Z}^n$7, extended to all $S^-\subseteq\mathbb{Z}^n$8 by range reduction:

$S^-\subseteq\mathbb{Z}^n$9

The optimization objectives are precision and cost, with real-valued max relative error

$S^\to$0

and float-valued max ULP error

$S^\to$1

The outer search uses precision $S^\to$2 or $S^\to$3 and cost as a $S^\to$4-D Pareto objective (Real et al., 2023).

Candidate programs are directed acyclic compute graphs with at most 100 vertices, consisting of one input vertex, zero or more coefficient vertices, internal binary-op vertices from $S^\to$5, and exactly one output vertex (Real et al., 2023). The explicit graph structure permits reuse of subexpressions and enables the discovery of non-polynomial patterns such as nested reciprocals and continued-fraction-style forms.

The search is bi-level. For a fixed graph structure, an inner CMA-ES loop optimizes coefficients over a training set of $S^\to$6 samples. An outer loop evolves a population of size $S^\to$7 using a multi-objective genetic-programming method called dNSGA-II, initialized with copies of the identity function $S^\to$8 (Real et al., 2023). Variation uses insertion, deletion, reconnection, or no-op mutations. Each child is then re-evaluated by coefficient optimization, JAX+XLA compilation, precision measurement on a hold-out set of size $S^\to$9 the training set, and throughput timing.

The reported empirical results are specific. For real $p\to q$0, a discovered 10-op approximation is stated to satisfy

$p\to q$1

which the paper characterizes as “essentially full double-precision accuracy with only 10 real operations” (Real et al., 2023). The comparison given is that “any 10-term Taylor or minimax polynomial is limited to $p\to q$2 error under the same cost.” For float32 $p\to q$3, the “top evolved candidate (12 float32 operations) achieves max ULP error $p\to q$4,” and is reported to be “about $p\to q$5 faster than the best Taylor/Padé/minimax baseline implemented in Horner form with FMAs” (Real et al., 2023).

The code-generation and compiler interaction are central to the float32 result. The evolved program is translated to JAX NumPy and JIT-compiled with XLA and LLVM. On “an Intel Skylake core,” the paper reports “$p\to q$6 billion evaluations/sec vs. $p\to q$7 billion for the best minimax polynomial—over $p\to q$8 faster—while staying under 1 ULP” (Real et al., 2023). The explanation given is that the discovered computation fuses into “two small fused computations,” which XLA classifies as “I/O bound,” whereas baseline implementations fuse into one larger computation classified as “CPU bound,” triggering a more expensive dispatch layer.

The paper also reports extension beyond the exponential. For log₂(x) on [1,2), “no classical baseline achieved \<3 ULP yet evolution finds \<1 ULP”; for erf(x) on [0,2] and a shifted Airy function on [0,1], evolved 5- to 20-op programs “match or exceed Padé/Chebyshev baselines by factors of 2–1 000 in precision for the same op count” (Real et al., 2023).

This version of AutoNumerics is therefore neither verification-oriented nor error-accounting-oriented. Its defining feature is black-box program discovery over a vast nonclassical design space, with the resulting approximants judged by precision and compilation behavior rather than by derivation from established approximation theory.

5. AutoNumerics as an autonomous PDE-solver pipeline

The most recent use of the name, "AutoNumerics: An Autonomous, PDE-Agnostic Multi-Agent Pipeline for Scientific Computing," frames AutoNumerics as a multi-agent system that converts natural-language PDE descriptions into transparent classical numerical solvers (Du et al., 19 Feb 2026). The framework consists of Formulator, Planner, Feature, Selector, Coder, Critic, and Reasoning agents, coordinated by a central dispatcher. The Formulator Agent parses the problem into a structured specification containing the PDE operator $p\to q$9, source $$\forall(p\to q)\in S^\to:\ (p\in S^+ \Rightarrow q\in S^+) \ \text{and}\ (q\in S^- \Rightarrow p\in S^-).$$0, boundary and initial conditions, and parameters. The Planner Agent enumerates candidate discretization and time-stepping schemes. The Feature Agent extracts numerical features such as linear versus non-linear, hyperbolic/parabolic/elliptic, and boundary types. The Selector Agent ranks candidate schemes by expected accuracy, stability margins, and computational cost. The Coder Agent generates executable solver code in Python/C++, the Critic Agent debugs it in two phases, and the Reasoning Agent interprets theoretical convergence or residual norms (Du et al., 19 Feb 2026).

The execution model is coarse-to-fine. In the coarse pre-check, the code is run on a downscaled mesh with $$\forall(p\to q)\in S^\to:\ (p\in S^+ \Rightarrow q\in S^+) \ \text{and}\ (q\in S^- \Rightarrow p\in S^-).$$1 to expose syntax errors, dimension mismatches, and shape errors. In the fine target run, the coarse solution $$\forall(p\to q)\in S^\to:\ (p\in S^+ \Rightarrow q\in S^+) \ \text{and}\ (q\in S^- \Rightarrow p\in S^-).$$2 is lifted to the fine grid by an interpolation operator

$$\forall(p\to q)\in S^\to:\ (p\in S^+ \Rightarrow q\in S^+) \ \text{and}\ (q\in S^- \Rightarrow p\in S^-).$$3

yielding

$$\forall(p\to q)\in S^\to:\ (p\in S^+ \Rightarrow q\in S^+) \ \text{and}\ (q\in S^- \Rightarrow p\in S^-).$$4

If the fine-grid run crashes or violates CFL/stability constraints, the Critic Agent adjusts $$\forall(p\to q)\in S^\to:\ (p\in S^+ \Rightarrow q\in S^+) \ \text{and}\ (q\in S^- \Rightarrow p\in S^-).$$5 or stencil parameters. A “Fresh Restart” regenerates code from scratch if either retry budget is exhausted (Du et al., 19 Feb 2026).

Verification is residual-based when no analytical solution is available. With discrete operator $$\forall(p\to q)\in S^\to:\ (p\in S^+ \Rightarrow q\in S^+) \ \text{and}\ (q\in S^- \Rightarrow p\in S^-).$$6 and source $$\forall(p\to q)\in S^\to:\ (p\in S^+ \Rightarrow q\in S^+) \ \text{and}\ (q\in S^- \Rightarrow p\in S^-).$$7, the pointwise residual is

$$\forall(p\to q)\in S^\to:\ (p\in S^+ \Rightarrow q\in S^+) \ \text{and}\ (q\in S^- \Rightarrow p\in S^-).$$8

and the relative residual norm is

$$\forall(p\to q)\in S^\to:\ (p\in S^+ \Rightarrow q\in S^+) \ \text{and}\ (q\in S^- \Rightarrow p\in S^-).$$9

with $A$0 to avoid division by zero. When an exact solution $A$1 is known, the relative error is

$A$2

For implicit relations $A$3, the framework defines

$A$4

Acceptance occurs when $A$5 or $A$6, with user-specified tolerance such as $A$7 (Du et al., 19 Feb 2026).

The scheme-selection rules reported in the paper map PDE structure to standard discretization families. The listed examples are: periodic hyperbolic $A$8 Fourier spectral + explicit RK4; periodic parabolic $A$9 Fourier spectral + IMEX; Dirichlet elliptic $+$0 Chebyshev spectral (Galerkin/Chebyshev-Tau); Dirichlet parabolic $+$1 finite difference or FEM with Crank-Nicolson or IMEX; and Dirichlet hyperbolic $+$2 high-order FD upwind or spectral (Du et al., 19 Feb 2026). This is significant because the framework explicitly positions itself against black-box neural solvers, emphasizing transparent solvers grounded in classical numerical analysis.

The experimental results are extensive. The paper reports experiments on 24 canonical and real-world PDE problems, with a CodePDE benchmark covering 1D Advection, 1D Burgers, 2D Reaction–Diffusion, 2D Compressible Navier–Stokes, and 2D Darcy Flow (Du et al., 19 Feb 2026). Under the reported settings, the Planner generates 10 plans, the top-5 are implemented, the retry limits are $+$3 and $+$4, fresh restarts are at most 2, the wall-clock cap per run is 120 s, and tolerances are $+$5 and $+$6.

For the CodePDE benchmark, the table lists the following geometric-mean values where available: FNO at $+$7, CodePDE (GPT-4.1) at $+$8, and AutoNumerics at $+$9 (Du et al., 19 Feb 2026). The paper states that AutoNumerics achieves “six orders of magnitude lower nRMSE than CodePDE and FNO.” On the selected 24-PDE table, representative entries include 1D Advection with error $\mathbb{R}$00 (L²) and runtime $\mathbb{R}$01 s, Poisson with error $\mathbb{R}$02 (L²) and runtime $\mathbb{R}$03 s, KdV with error $\mathbb{R}$04 (L²) and runtime $\mathbb{R}$05 s, Gray–Scott without analytical solution with residual $\mathbb{R}$06 and runtime $\mathbb{R}$07 s, Biharmonic with error $\mathbb{R}$08 (L²), and 5D Helmholtz with error $\mathbb{R}$09 (L²) (Du et al., 19 Feb 2026). The paper further reports that convergence orders match “2nd order for FD” and “spectral exponential rate for smooth periodic problems,” with runtimes from 18 s to 230 s.

The limitations are explicit: the framework “struggles with very high-order (4th order) or high-dimensional ($\mathbb{R}$10D) PDEs,” is “restricted to regular Cartesian domains,” is “dependent on a single LLM (GPT-4.1),” and “lacks formal proofs of convergence or stability for generated code” (Du et al., 19 Feb 2026). These caveats are integral to understanding the scope of the reported autonomy.

6. Antecedents and adjacent formulations

Monniaux’s "Automatic modular abstractions for template numerical constraints" provides a notable antecedent to later automation efforts in numerical reasoning (0909.4013). The method starts from a template linear-constraint domain defined by fixed linear forms

$\mathbb{R}$11

and template formula

$\mathbb{R}$12

Its concretization is

$\mathbb{R}$13

and the abstract counterpart is

$\mathbb{R}$14

The pair $\mathbb{R}$15 forms a Galois connection (0909.4013).

For loop-free code with semantics $\mathbb{R}$16, soundness of a postcondition template $\mathbb{R}$17 requires

$\mathbb{R}$18

Quantifier elimination in Linear Real Arithmetic is then used to derive, for each output parameter, the least sound bound. For single loops, the invariant template $\mathbb{R}$19 is constrained by

$\mathbb{R}$20

again with least inductive invariants extracted by quantifier elimination (0909.4013). The output is compiled into piecewise-linear if-then-else code computing optimal abstract transformers or loop invariants as functions of precondition parameters.

The toolchain parses a small C-like language into logical formulas, constructs soundness or invariant conditions, invokes quantifier-elimination engines, and pretty-prints the resulting C code (0909.4013). The paper reports that generating a transformer for a 20-line numeric library block takes “1–2 s over reals” and “10–20 s when encoding IEEE-754 rounding,” and notes that commercial packages often run out of memory on the same tasks. Its limitations include the need for a global view of the whole loop or function, the doubly-exponential worst-case cost of quantifier elimination, and the fact that template linear forms must be chosen a priori.

A different adjacent formulation comes from semiautomatic rings and grids (Gao et al., 2021). There, an “automatic numeric system” consists of a countable dense subring of $\mathbb{R}$21 with automatic addition and comparisons and semiautomatic multiplication by constants. The canonical example is the $\mathbb{R}$22-adic rationals

$\mathbb{R}$23

encoded as finite signed base-$\mathbb{R}$24 strings (Gao et al., 2021). The paper then constructs richer dense semiautomatic grids containing algebraic numbers such as $\mathbb{R}$25 or $\mathbb{R}$26, enabling automatic representation of geometric transformations requiring those constants.

The geometric consequences are precise. Translation by a fixed vector is automatic. Rotation by $\mathbb{R}$27 requires

$\mathbb{R}$28

and is automatic in grids containing $\mathbb{R}$29 and $\mathbb{R}$30; rotation by $\mathbb{R}$31 requires $\mathbb{R}$32 and $\mathbb{R}$33 (Gao et al., 2021). The paper also states that in every semiautomatic grid, the collection of all affine lines is not an automatic family, while each individual triangle with corner points in the grid is regular. As a corollary, nondegenerate equilateral triangles cannot have all three vertices in pure $\mathbb{R}$34-adic rationals $\mathbb{R}$35, but can be realized in the constructed $\mathbb{R}$36-grid.

These antecedents are not labeled AutoNumerics in their titles, but they clarify that automated numerical reasoning has long included both symbolic synthesis and automata-theoretic representation. This suggests that the later uses of the name occupy different points in an already heterogeneous design space.

7. Comparative interpretation, misconceptions, and limitations

A common misconception would be to treat AutoNumerics as a single method with one objective function, one implementation strategy, or one research community. The literature does not support that interpretation. The invariant-synthesis system uses ICE samples, abstract domains, and SMT checks (Bouajjani et al., 2022); the precision-estimation framework propagates exact_bits counters through IEEE-754 operations (Netay, 2024); AutoNumerics-Zero performs evolutionary search over compute graphs compiled with JAX+XLA (Real et al., 2023); and the PDE pipeline orchestrates specialized LLM agents to construct classical solvers from natural language (Du et al., 19 Feb 2026).

A second misconception would be that automation in numerics is synonymous with neural operators or black-box learned surrogates. The PDE-oriented AutoNumerics explicitly contrasts itself with “black-box neural solvers” and instead generates finite difference, finite element, spectral, and related classical schemes (Du et al., 19 Feb 2026). Conversely, AutoNumerics-Zero shows that automated search can still be black-box at the design stage even when the discovered output is an explicit analytical compute graph (Real et al., 2023). The contrast is therefore between how the method is found and what form the resulting numerical object takes.

The limitations across the variants are also structurally different. Numerical invariant synthesis depends on the chosen abstract domain and on the ability of separator-derived constraints to support a sufficient attribute pool (Bouajjani et al., 2022). Automatic precision estimation is conservative by design and introduces measurable storage and runtime overheads (Netay, 2024). AutoNumerics-Zero relies on large-scale evolutionary search and compiler-specific performance behavior, and its strongest reported speedups exploit XLA/LLVM compilation paths rather than solely approximation-theoretic superiority (Real et al., 2023). The PDE pipeline lacks formal proofs of convergence or stability for generated code and is restricted to regular Cartesian domains (Du et al., 19 Feb 2026). Monniaux’s synthesis method, meanwhile, is constrained by quantifier-elimination cost and by a priori template selection (0909.4013).

A plausible implication is that the most stable unifier of the AutoNumerics label is not a specific algorithmic core, but a methodological stance: numerical structures that were once hand-designed are increasingly being generated, selected, or audited automatically under explicit semantic criteria. Those criteria differ markedly—inductiveness, guaranteed correct bits, ULP error, nRMSE, residual norms, or automaton recognizability—but each framework attempts to turn expert numerical craftsmanship into a reproducible computational pipeline.