The VHDL standard's float package gets sqrt(2) wrong

A hardware arithmetic audit uncovered precision defects in the official floating-point package for the IEEE 1076 VHDL standard, where basic operations return incorrect values under default settings. Detailed in github.com's report, evaluating the square root of 2.0 in double precision produces a result one unit in the last place too low. The error affects IEEE Std 1076-2008 and 1076-2019, whose standard packages ship unmodified in open-source simulators like GHDL and NVC.

This analysis covers the root causes of the calculation errors, the failure rates discovered across standard numeric formats, and the corrected package patches submitted to language maintainers. It is written for FPGA designers, verification engineers, and toolchain developers who rely on VHDL simulation libraries for hardware modeling. Readers will find exact reproduction steps, failing operand examples, and verification metrics comparing standard functions against exact rational arithmetic.

01Root iterations and missing sticky bits

The core square root routine in the VHDL standard library relies on a fixed count of Newton steps rather than an exact bit-by-bit extraction. Because the implementation lacks a sticky bit to track discarded fractional remainder data, the calculation truncates precision before final rounding. In contrast, taking the root one bit at a time using one subtraction per bit yields an exact result.

Arithmetic errors extend beyond square roots into multiply-accumulate operations and conversions. In test evaluations, evaluating mac(-6, -1, -5) produces -15 instead of the correct value 1. Converting a double-precision 2^-512 to single precision unexpectedly returns 1.0, while taking the reciprocal of 2^127 in single precision returns positive infinity.

Testing arithmetic against exact fractional references reveals systematic rounding failures across thousands of standard calls. In an exhaustive check of a 7-bit format with three exponent and three fraction bits, 1,271,619 calls failed to round correctly out of 8,692,326 total runs. Furthermore, the related fixed-point package division routine fails on 41,026 out of 751,418 test calls.

02Standard packages across GHDL and NVC

The flawed code resides directly inside ieee/float_generic_pkg-body.vhdl, an IEEE-maintained source file spanning 5,712 lines under revision 1220. This canonical package is distributed on opensource.ieee.org/vasg/Packages and the IEEE-P1076 working group repository. Open-source simulation engines incorporate these identical files without downstream modifications.

Validation runs confirmed identical failure patterns across independent toolchains, specifically GHDL version 7.0.0-dev and NVC version 1.23.1. When fed identical test vectors, both simulators generated matching output files containing identical numeric errors. The audit evaluated simulator behavior against exact.py, a reference script that performs fractional arithmetic before applying target rounding modes.

Simulators can reproduce the issues using standard commands without specialized test harnesses. Compiling plain.vhd with standard 2008 compliance flags allows GHDL to execute the plain test suite directly from the terminal. The execution prints explicit failure indicators for mismatched values before halting with a boundary error inside the package body.

03Plain testbench calls and simulator crashes

Under default package settings with round-nearest, three guard bits, and denormal support enabled, standard calls return corrupted values. Evaluating the square root of 14.0 produces hex 406F7750 instead of 406F7751, while the square root of 2.0^-149 yields 4.2e-22 rather than 3.7e-23. In addition, mac(1.0, 1.0, -infinity) produces positive infinity instead of negative infinity.

Format resize operations exhibit extreme edge-case instability across multiple exponent bands. Converting double-precision 2.0^-512 to float32 produces 1.0 instead of 0.0, while 2.0^600 converts to 2.0^88 instead of infinity. Attempting to resize 2.0^-650 halts the simulation entirely due to an index out of bounds error at line 358 of the package body.

Failures appear consistently across ordinary integer inputs rather than isolated mathematical anomalies. Evaluating square roots across integers 1 to 4,096 returns an incorrect value for 256 inputs in float32 and 264 inputs in float64. Out of 4,913 simple multiply-accumulate triples, 290 float32 calculations return invalid values under standard default rounding.

04Millions of evaluations across precision formats

Stress testing applied 10,000 calls per operation across four standard rounding modes targeting arithmetic boundaries like denormals, exponent limits, exact ties, and signed zeros. In float32 testing, the package logged 81,068 incorrect results and 272 simulator terminations across 859,771 evaluated calls. Testing float64 produced 86,081 wrong values and 445 simulation crashes across 857,234 calls.

Extended formats showed comparable levels of calculation corruption under structured test conditions. Evaluations in float16 registered 113,980 incorrect returns across 866,782 calls, while float128 logged 90,744 calculation errors alongside 585 complete simulator terminations. Dedicated scaling and division routines like scalb and dividebyp2 repeatedly delivered corrupted outputs across thousands of trial vectors.

Exhaustive evaluation of reduced-width formats confirmed deep algorithmic flaws in rounding logic. Across a 10-bit format with four exponent and five fraction bits, 883,716 operations out of 17,340,682 calls were incorrect. Single-operand float16 functions produced 266,342 errors out of 2,998,266 tests, while round-nearest mac failed on 243,520 of 2,097,152 possible triples.

05Working group patches and passing testbenches

Analysis of error distances under round-nearest modes demonstrates severe deviations from mathematical precision. While correctly rounded operations must remain within 0.5 units in the last place, over 1,900 float32 calls deviated by more than 1,000 units. More than 2,000 float32 calls returned incorrect signs, and 1,897 delivered infinity or NaN instead of finite values.

The audit produced a comprehensive correction that replaces flawed iterative routines with bit-exact operations. Across four extensive validation test sets encompassing 32,470,838 total calls, every single test vector passed successfully once the patch was applied. Both GHDL and NVC compiled the patched source and generated identical, fully verified output files.

Formal defect reports have been submitted to upstream tool maintainers and standards bodies to coordinate vendor updates. The issue is tracked under IEEE P1076 working group ticket VHDL-Issues 388, GHDL issue 3486, NVC issue 1758, and fphdl issue 4. The distribution repository includes complete source patches alongside testbenches to allow teams to verify local fixes immediately.