# Verifying numerical code

A **numerical bug** is a program that computes the wrong number and returns it
without failing. There is no crash and no exception, only a plausible wrong
value, so the difficulty is finding an authority to check the answer against
rather than noticing that something went wrong.

Di Franco et al. categorize numerical bugs in released libraries into four
kinds (Di Franco et al. 2017)[^difranco2017]. The kinds need different methods:

- **Correctness** — the mathematics itself is wrong: a formula transcribed
  incorrectly, a mismatched array dimension, an assumption the compiler was
  free to break.
- **Special value** — a NaN, an infinity, a signed zero, or a subnormal enters
  and propagates instead of being handled, often out of a domain error such as
  the logarithm of zero.
- **Convergence** — an iterative method stops at the wrong point, oscillates,
  or does not terminate.
- **Accuracy** — rounding or truncation loses enough precision to change the
  answer.

Underneath all four sits the choice of representation, which removes failure
modes instead of detecting them.

## Choosing the representation

**A fixed-width integer is the integers modulo 2ⁿ.** Wraparound, truncating
conversion, mixing signed with unsigned, and shifting past the bit width all
yield a defined-but-wrong value. A truncating conversion of a floating-point
value to a 16-bit integer is what destroyed Ariane 5 flight 501 (Dietz et al. 2015)[^dietz2015].
C and C++ are the exception that matters: signed overflow there is undefined
rather than wrapping, so an optimizer is entitled to assume it cannot happen
and to delete the check written to catch it. The escape is an
arbitrary-precision integer (Python's `int`, Java's `BigInteger`,
[GMP](https://gmplib.org/)) where the cost is acceptable, and explicit fixed-width
arithmetic where it is not: Rust's `checked_*`, `wrapping_*`, and
`saturating_*` families with overflow panics in debug builds, Swift's trapping
operators alongside `&+`, Zig's separation of `+` from `+%`. What those buy is
not the absence of overflow but the absence of *unintended, silent* overflow.

**A binary float is a finite subset of the reals, rounded at every step.**
Addition is not associative, so a sum depends on the order of its terms, and
subtracting nearly equal quantities cancels the leading digits and promotes
rounding noise into the result.

**A binary float is not a decimal.** A fraction is exact in base *b* only when
every prime factor of its denominator divides *b*. Ten factors into 2 and 5, so
a tenth is exact in decimal; two factors into itself alone, so a tenth repeats
in binary and is stored as the nearest float instead. Adding `0.1` and `0.2` and
printing `0.30000000000000004` is the standard demonstration, and it reproduces
in most languages (Wiffin 2012)[^wiffin2012]. Money is the usual casualty, and the escape is
a different representation rather than more precision: a decimal type (Java's
`BigDecimal`, Python's `decimal`, SQL `NUMERIC`) or storage in integer minor
units.

**A number carries no unit.** Meters and feet share a type, so a dimensional
mistake is invisible to the compiler. F# has units of measure in the language;
elsewhere a newtype per quantity reaches the same place.

## Getting the mathematics right

A formula mis-transcribed from a paper is wrong for every input, which would
make it easy to catch if anything existed to compare against, and frequently
nothing does: the program *is* the definition of the value it computes. Four
ways to manufacture an authority:

- **An independent implementation.** [Differential
  testing](https://quality.stereobooster.com/differential-testing.md) against another library, a
  computer-algebra system, or an exact symbolic computation over small inputs.
  A deliberately naive implementation, too slow to ship, is a legitimate
  reference for the cases it can handle.
- **A relation instead of a value.** Where no reference exists, a [metamorphic
  relation](https://quality.stereobooster.com/metamorphic-testing.md) constrains the output without
  naming it: scale invariance, linearity, symmetry, a forward transform undone
  by its inverse, an identity that must hold between two calls. The relations
  metamorphic testing was introduced with are numerical ones.
- **A residual.** For anything that solves an equation, substituting the answer
  back in is an oracle that is both cheaper than the solve and independent of
  it. This checks the answer rather than the algorithm, which is what makes it
  survive a rewrite.
- **Reading.** A transcription error against a published formula is found by
  comparing source to paper, which is [code
  review](https://quality.stereobooster.com/code-review.md) by someone who knows the mathematics, not
  a test.

## Handling the special values

Special values are the second-largest category and the cheapest to reach,
because the inputs that produce them are a short, known list.
[Property-based testing](https://quality.stereobooster.com/property-based-testing.md) supplies them by
default: [Hypothesis](https://hypothesis.readthedocs.io/)'s `floats()` emits NaN, both infinities,
both zeros, and subnormals unless told otherwise, as do
[fast-check](https://fast-check.dev/) and [jqwik](https://jqwik.net/). A generator narrowed to
plausible-looking values excludes exactly the inputs this category is made
of.

IEEE 754 specifies floating-point exception conditions and their default
handling (IEEE 2019)[^ieee754_2019]; where a platform allows those conditions to trap
instead (`feenableexcept` in glibc, on invalid, overflow, and divide-by-zero),
a NaN that would otherwise propagate silently becomes a fault at the
instruction that created it, where diagnosis is cheap. A
[precondition](https://quality.stereobooster.com/contracts-and-runtime-assertions.md) stating the
mathematical domain of a function (a non-negative argument to a square root, a
positive one to a logarithm) turns a silent NaN into a reported contract
violation at the boundary. And [abstract
interpretation](https://quality.stereobooster.com/abstract-interpretation.md) over a domain that tracks
these values finds them statically rather than by sampling.

NaN compares unequal to everything including itself, so it violates the
contract every comparison-based sort and every ordered container assumes. A
collection holding one can produce a corrupted order or an infinite loop far
from the arithmetic that created it.

## Convergence and termination

An iterative method has two failure modes that a returned value cannot
distinguish: it stopped early, or it stopped at the wrong place. Neither shows
up as an exception unless the code raises one, so the first requirement is that
every iteration carry a cap and report exhaustion as a failure rather than
returning its last iterate. [Termination
analysis](https://quality.stereobooster.com/termination-analysis.md) proves the stronger property that
the loop ends at all, which matters where the iteration count is data-dependent
and unbounded in principle.

The oracle for whether it stopped in the right place is the residual, not the
iterate: an iterate that stops moving indicates the method converged to
*something*. Where the input governs whether convergence is even possible, the
condition number is the quantity that says so, and estimating it is a
precondition on the input rather than a check on the output.

## Bounding the accuracy

Rounding error is continuous, so the question is not whether it occurred but
how large it can grow over the input range, and the answer is a bound.
[Abstract interpretation](https://quality.stereobooster.com/abstract-interpretation.md) supplies one
from numeric domains built for the purpose: intervals for a coarse range, and
affine or polyhedral domains where correlations between variables matter and an
interval bound would explode. [IKOS](https://github.com/NASA-SW-VnV/ikos), [MOPSA](https://mopsa.lip6.fr/), Frama-C's
EVA plugin ([Frama-C](https://www.frama-c.com/)), and the commercial Astrée run these.

A second family targets round-off specifically, reporting a bound on the
difference between the floating-point result and the real-valued one:
[Gappa](https://gappa.gitlabpages.inria.fr/), [FPTaylor](https://github.com/soarlab/FPTaylor), [Daisy](https://gitlab.mpi-sws.org/AVA/daisy),
[PRECiSA](https://github.com/nasa/PRECiSA), and the commercial Fluctuat, with
[FPBench](https://fpbench.org/) as the shared benchmark and exchange format. Several
emit a certificate a proof assistant can check, which moves the result from a
tool's assertion to [theorem proving](https://quality.stereobooster.com/theorem-proving.md): Gappa
produces [Coq](https://rocq-prover.org/) proofs, FPTaylor produces HOL Light ones.

Sampling reaches what a bound cannot. Because the error is a continuous
function of the input, finding the input that maximizes it is an optimization
problem, which is [search-based software
testing](https://quality.stereobooster.com/search-based-software-testing.md) with a numeric objective;
[Herbie](https://herbie.uwplse.org/) works this way, sampling to localize error in an
expression and then searching a rule database for a more accurate equivalent
form. [Verificarlo](https://github.com/verificarlo/verificarlo) and [CADNA](http://cadna.lip6.fr/) take the opposite
approach and perturb rounding across repeated runs, reporting how many digits
of the output stay fixed, which estimates the significance of a result with no
reference implementation at all.

Where a reference is available, it is the sharpest oracle of the group: the
same computation carried out in arbitrary precision ([MPFR](https://www.mpfr.org/),
`mpmath`, Python's `fractions`) gives the correctly-rounded answer, compared in
units in the last place rather than as an absolute difference. An epsilon
chosen because it made the test pass is not a tolerance; the number has to come
from an error analysis.

## Evidence

Di Franco et al. inspected 828 bug reports across NumPy, SciPy, LAPACK, the GNU
Scientific Library, and Elemental, identified 269 as numerical, and sorted
them: correctness 99, special value 76, convergence 56, accuracy 38
(Di Franco et al. 2017)[^difranco2017]. Correctness bugs, the largest group, are transcription and
dimension errors rather than artifacts of floating point, so a toolset aimed
only at rounding addresses a seventh of the measured population.

For the integer side, Dietz et al. built the IOC dynamic checker and measured
overflow in real C and C++: more than 200 distinct locations in the SPEC
CINT2000 benchmarks, most of them deliberate wraparound rather than defects,
plus undefined signed overflows found and reported in SQLite, PostgreSQL,
SafeInt, GNU MPC and GMP, Firefox, LLVM, Python, BIND, and OpenSSL
(Dietz et al. 2015)[^dietz2015]. Because most of the measured overflows are intentional, a checker
that cannot separate deliberate wraparound from accidental overflow reports
both.

Three results measure tools rather than bugs. Herbie improved accuracy on every
example drawn from a numerical-methods textbook, by up to 60 bits on some, at a
median 40% runtime overhead (Panchekha et al. 2015)[^panchekha2015]. Solovyev et al. ran six
round-off bound tool families over 24 benchmark examples on one machine with
each tool's options tuned; the authors, who wrote one of the six, report their
own tool computing the tightest bound on 23 of the 24, against 3 of 20 for the
next best (Solovyev et al. 2018)[^solovyev2018]. DEBAR, an abstract
interpreter for special-value bugs in neural-network architectures, detected
every bug in a known-bug dataset with no false positives, and on real-world
architectures reported 529 warnings of which 299 were true positives
(Zhang et al. 2020)[^zhang2020].

## Referenced by

- [Metamorphic testing](https://quality.stereobooster.com/metamorphic-testing.md) · Methods

## References

[^difranco2017]: Di Franco, Anthony, Hui Guo, and Cindy Rubio-González. 2017. "[A Comprehensive Study of Real-World Numerical Bug Characteristics](https://web.cs.ucdavis.edu/~rubio/includes/ase17.pdf)." *Proceedings of the 32nd IEEE/ACM International Conference on Automated Software Engineering (ASE 2017)*, 509–19. <https://doi.org/10.1109/ASE.2017.8115662>.
[^dietz2015]: Dietz, Will, Peng Li, John Regehr, and Vikram Adve. 2015. "[Understanding Integer Overflow in C/C++](https://www.cs.utah.edu/~regehr/papers/tosem15.pdf)." *ACM Transactions on Software Engineering and Methodology* 25 (1): 1–29. <https://doi.org/10.1145/2743019>.
[^wiffin2012]: Wiffin, Erik. 2012. *[Floating Point Math](https://0.30000000000000004.com/)*. 0.30000000000000004.com. <https://0.30000000000000004.com/>.
[^ieee754_2019]: IEEE. 2019. *[IEEE Std 754-2019 — IEEE Standard for Floating-Point Arithmetic](https://standards.ieee.org/ieee/754/6210/)*. IEEE. <https://doi.org/10.1109/IEEESTD.2019.8766229>.
[^panchekha2015]: Panchekha, Pavel, Alex Sanchez-Stern, James R. Wilcox, and Zachary Tatlock. 2015. "[Automatically Improving Accuracy for Floating Point Expressions](https://herbie.uwplse.org/pldi15-paper.pdf)." *Proceedings of the 36th ACM SIGPLAN Conference on Programming Language Design and Implementation (PLDI 2015)*, 1–11. <https://doi.org/10.1145/2737924.2737959>.
[^solovyev2018]: Solovyev, Alexey, Marek S. Baranowski, Ian Briggs, Charles Jacobsen, Zvonimir Rakamarić, and Ganesh Gopalakrishnan. 2018. "[Rigorous Estimation of Floating-Point Round-Off Errors with Symbolic Taylor Expansions](https://soarlab.org/papers/2018_toplas_sbbjrg.pdf)." *ACM Transactions on Programming Languages and Systems* 41 (1): 1–39. <https://doi.org/10.1145/3230733>.
[^zhang2020]: Zhang, Yuhao, Luyao Ren, Liqian Chen, Yingfei Xiong, Shing-Chi Cheung, and Tao Xie. 2020. "[Detecting Numerical Bugs in Neural Network Architectures](https://dl.acm.org/doi/pdf/10.1145/3368089.3409720)." *Proceedings of the 28th ACM Joint European Software Engineering Conference and Symposium on the Foundations of Software Engineering (ESEC/FSE 2020)*, 826–37. <https://doi.org/10.1145/3368089.3409720>.

## Acronyms

- HOL — higher-order logic
- IOC — integer overflow checker