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)¹. 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)². 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) 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)³. 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 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 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 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 supplies them by default: Hypothesis's floats() emits NaN, both infinities, both zeros, and subnormals unless told otherwise, as do fast-check and jqwik. 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)⁴; 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 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 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 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 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, MOPSA, Frama-C's EVA plugin (Frama-C), 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, FPTaylor, Daisy, PRECiSA, and the commercial Fluctuat, with FPBench 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: Gappa produces Coq 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 with a numeric objective; Herbie works this way, sampling to localize error in an expression and then searching a rule database for a more accurate equivalent form. Verificarlo and CADNA 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, 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)¹. 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)². 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)⁵. 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)⁶. 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)⁷.

Referenced by

• Metamorphic testing · Methods

References

---

1. Di Franco, Anthony, Hui Guo, and Cindy Rubio-González. 2017. "A Comprehensive Study of Real-World Numerical Bug Characteristics." Proceedings of the 32nd IEEE/ACM International Conference on Automated Software Engineering (ASE 2017), 509–19. https://doi.org/10.1109/ASE.2017.8115662. ↩↩

2. Dietz, Will, Peng Li, John Regehr, and Vikram Adve. 2015. "Understanding Integer Overflow in C/C++." ACM Transactions on Software Engineering and Methodology 25 (1): 1–29. https://doi.org/10.1145/2743019. ↩↩

3. Wiffin, Erik. 2012. Floating Point Math. 0.30000000000000004.com. https://0.30000000000000004.com/. ↩

4. IEEE. 2019. IEEE Std 754-2019 — IEEE Standard for Floating-Point Arithmetic. IEEE. https://doi.org/10.1109/IEEESTD.2019.8766229. ↩

5. Panchekha, Pavel, Alex Sanchez-Stern, James R. Wilcox, and Zachary Tatlock. 2015. "Automatically Improving Accuracy for Floating Point Expressions." Proceedings of the 36th ACM SIGPLAN Conference on Programming Language Design and Implementation (PLDI 2015), 1–11. https://doi.org/10.1145/2737924.2737959. ↩

6. 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." ACM Transactions on Programming Languages and Systems 41 (1): 1–39. https://doi.org/10.1145/3230733. ↩

7. Zhang, Yuhao, Luyao Ren, Liqian Chen, Yingfei Xiong, Shing-Chi Cheung, and Tao Xie. 2020. "Detecting Numerical Bugs in Neural Network Architectures." 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. ↩