Audit of VHDL ieee.float_pkg: 1,271,619 of 8,692,326 calls in a 7-bit format not correctly rounded
An exact-arithmetic check found rounding errors and simulator stops in the IEEE floating-point package that GHDL and NVC ship unchanged; the authors published a patch.
The authors of an audit published on github.com report that ieee.float_pkg, the floating-point package of the VHDL standard IEEE Std 1076-2008 and 1076-2019, returns results that are not correctly rounded. IEEE publishes the file, and GHDL and NVC ship it unchanged. For all operands and pairs in a 7-bit format (3 exponent bits, 3 fraction bits), plus all multiply-accumulate triples, the authors ran 8,692,326 calls, and 1,271,619 were not correctly rounded github.com
github.comgithub.comOpen source ↗.
The simplest case is the square root: sqrt(2.0) in double precision returns a value one unit in the last place too low, by default. It lies 0.565 of a unit from the root, while the correctly rounded value is 0.435 from it. Multiply-accumulate (-6)(-1)+(-5) returns -15 instead of 1; converting 2^-512 from double to single returns 1.0; the reciprocal of 2^127 in single returns infinity.
How the check was done
The check treats operands as exact fractions, performs each operation exactly and rounds once before comparing. The runs were made in GHDL 7.0.0-dev and NVC 1.23.1, and on the same vectors both wrote identical result files. Each operation, format and rounding style got 10,000 calls on operands chosen for places where arithmetic breaks: denormals, both ends of the exponent range, exact ties, signed zeros, infinities.
In float32 this gave 81,068 wrong results and 272 simulator stops in 859,771 calls; in float64, 86,081 wrong and 445 stops in 857,234 github.com
github.comgithub.comOpen source ↗. The ordinary operands are affected too: of the integers 1 to 4,096, sqrt is one unit off for 256 in float32 and for 264 in float64. The minimal test file has 25 calls: 24 print WRONG, and the 25th stops the simulator with an index out of bounds at line 358 of the package body.
Causes and patch
The authors trace the errors to specific lines of the IEEE file. In normalize (line 422) the exponent is cut to two bits more than the target format, so from float64 to float32 it is taken modulo 1,024. In divide (line 1243) and reciprocal (line 1099) the sticky bit passed to rounding is the constant '1', so exact ties go to the wrong neighbour. In reciprocal the result exponent does not fit in the top binade, and in float32 every argument from 2^127 up returns infinity. The square root uses a fixed count of Newton steps and no sticky bit.
The repository holds the exact reference, the testbench, the causes with line numbers and a patch. With the patch, all 32,470,838 calls from four test sets pass. The authors also report that division in fixed_pkg gets 41,026 of 751,418 results wrong and give a patch for it. The issues were reported to the IEEE P1076 working group (VHDL-Issues 388), ghdl 3486, nvc 1758 and fphdl 4.
