IEEE's VHDL float package mis-rounds sqrt(2), crashes simulators; GHDL and NVC affected

The floating-point package that ships with the VHDL standard computes sqrt(2.0) one unit in the last place below the correctly rounded result in double precision, according to an audit published in a GitHub repository and shared on Hacker News's front page. Because the defective file is published by IEEE itself as part of IEEE Std 1076-2008 and 1076-2019, the error reaches users through GHDL and NVC, both of which distribute the package unchanged.

The headline defect

Under the package's own defaults — round-to-nearest, three guard bits, denormal support on — the audit reports that sqrt(2.0) in float64 produces the hex value 3FF6A09E667F3BCC, while the exactly rounded answer is 3FF6A09E667F3BCD. The float32 sqrt(14.0) is likewise one bit low. The audit attributes the problem to an implementation that refines an initial root with a fixed number of Newton steps and keeps no sticky bit, and notes that a digit-by-digit method — one subtraction per result bit, the same procedure as longhand square root on paper — is exact by construction.

It also takes nothing exotic to trip the bug. Of the integers 1 through 4,096, the audit counts 256 whose float32 square root comes back one unit off and 264 whose float64 result does. Among 4,913 small multiply-add triples, 290 return the wrong float32.

How the results were checked

The audit's reference implementation treats every operand as a sign plus an exact rational fraction, performs each operation exactly, and rounds once into the destination format in whichever rounding style the call names. Simulator output is compared against that reference.

On a deliberately tiny 7-bit format (three exponent bits, three fraction bits), every operand, operand pair and multiply-add triple was enumerated: 8,692,326 calls, of which 1,271,619 were not correctly rounded. Larger formats were probed with 10,000 targeted calls per operation, format and rounding mode, aimed at the places arithmetic breaks — denormals, both ends of the exponent range, values that round into the next binade, exact ties, signed zeros and infinities. Even under round-to-nearest, mac misses 2,671 of 10,000 float32 calls and scalb misses 2,653. Directed rounding modes fare far worse: float32 divide in round_inf is wrong 2,395 times out of 10,000, and multiply in round_zero 1,102 times. Add, subtract and multiply are exactly rounded under round-to-nearest, with zero misses in 10,000, yet they too fail in the directed modes.

Errors beyond the last bit

Several results are grossly wrong rather than off by a single bit. mac(-6.0, -1.0, -5.0) returns -15.0 where the exactly rounded result is 1.0, and mac(1.0, 1.0, -infinity) returns +infinity instead of -infinity. Shrinking 2.0**-512 from float64 to float32 yields 1.0 instead of 0.0, and the reciprocal of 2.0127 in float32 returns +infinity instead of 2.0-127. Conversions drift as well: to_signed(0.50000006, 32) gives 0 instead of 1, and to_unsigned(199.78125, 8) gives 199 instead of 200.

Worst, some inputs halt the tool entirely. resize from float64 to float32 returns +infinity for every power of two from 2^-1074 to 2^-1024, returns values with the exponent shifted by 512 for other ranges, and stops the simulator with an index-out-of-bounds error at line 358 of the package body for inputs between 2^-662 and 2^-639. The audit's plain.vhd testbench needs only GHDL: after compiling and running with --std=08, 24 lines print WRONG and the 25th call stops the simulation.

Patches and upstream reports

The audited source is ieee/float_generic_pkg-body.vhdl, 5,712 lines at revision 1220, from the release branch of the IEEE's Packages repository. The working group's develop branch carries identical bytes, and the audit's runs on GHDL 7.0.0-dev and NVC 1.23.1 wrote identical result files wherever both ran the same vectors — agreement between tools, but agreement on the wrong answers. The audit publishes a patch, after which all 32,470,838 calls across four test sets are correctly rounded. The companion fixed-point package does not escape either: fixed_pkg division is wrong in 41,026 of 751,418 calls, with its own patch. The findings were reported to the IEEE P1076 working group as VHDL-Issues 388, and filed against the simulators as GHDL issue 3486 and NVC issue 1758.

Why it matters

A defect in a standard's own package is worse than a bug in one tool: anything that ships the IEEE-published file inherits it, and engineers cross-checking GHDL against NVC would see matching output and reasonably conclude it was correct. The failures go beyond last-bit pedantry — there is a sign error in multiply-add, conversions off by whole values, and inputs that crash the simulator — so testbenches and arithmetic models built on float_pkg may have been silently validating against wrong references. The audit also shows the method's value: exact rational arithmetic plus exhaustive testing on a tiny 7-bit format surfaced defects that years of shipping did not, and demonstrated that correctly rounded results, such as digit-by-digit square root, are achievable in a standard library.