We established sqrt(2.0) in the standard's float_pkg returns one component in the final location too low in twice precision by default. This is the IEEE 1076-2008 and 1076-2019 VHDL floating-point package; IEEE publishes the file, GHDL and NVC container it unchanged.
The base uses a fixed figure of Newton steps and no sticky bit. But taking it one bit at a period from the definition is exact, alongside a subtraction per bit, the identical item as quadrate base by hand on paper.
Our inspect treats operands as exact fractions, does all procedure exactly and rounds formerly before comparing. mac(-6, -1, -5), or (-6)(-1)+(-5), returns -15 alternatively of 1. Converting 2^-512 twice to sole returns 1.0; reciprocal of 2^127 in sole returns infinity. For all operands and pairs in a 7-bit format (3 exponent bits, 3 fraction bits), affirmative all mac triples, we ran 8,692,326 calls, and 1,271,619 weren't correctly rounded.
The exact reference, testbench, causes alongside origin row numbers and place are in the folder. All 32,470,838 calls from four test sets continue alongside patch; GHDL and NVC create identical outcome records on average sets. fixed_pkg division gets 41,026 of 751,418 incorrect and has its place there. Reported to IEEE P1076 operating collection (VHDL-Issues 388), ghdl 3486, nvc 1758.
To see it yourself, alongside GHDL and nothing else:
ghdl -a --std=08 plain.vhd ghdl --elab-run --std=08 plain
The remainder of this leaf is the record: all incorrect call, the counts, the causes by line, the patch.
ieee.float_pkg is the floating-point bundle of the VHDL standard, IEEE Std 1076-2008 and 1076-2019. The origin is the IEEE's own: ieee/float_generic_pkg-body.vhdl on the publish branch of opensource.ieee.org/vasg/Packages, 5,712 lines, $Revision: 1220 $. The operating group's create branch at gitlab.com/IEEE-P1076/Packages has the identical bytes, and GHDL and NVC container them unchanged. The runs current are in GHDL 7.0.0-dev and NVC 1.23.1, and anywhere the two ran the identical vectors they wrote the identical outcome file.
The bundle is set beside exact arithmetic. An operand is a sign and a fraction of two integers, the procedure is done on the fractions, and the outcome is rounded formerly into the format in the rounding manner the call names. exact.py is that reference.
With the package's own settings: round_nearest, 3 defender bits, denormal numbers on.
call returns exactly rounded sqrt (2.0) in float64 3FF6A09E667F3BCC 3FF6A09E667F3BCD sqrt (14.0) 406F7750 406F7751 sqrt (2.0**-149) 4.2e-22 3.7e-23 mac (-6.0, -1.0, -5.0) -15.0 1.0 mac (-5.0, 1.0, 4.0) 15.0 -1.0 mac (1.0, 1.0, -infinity) +infinity -infinity mac (1 + 2.0**-23, 1 - 2.0**-23, -1.0) -2.0**-26 -2.0**-46 resize (2.0**-512 in float64, 8, 23) 1.0 0.0 resize (2.0**600 in float64, 8, 23) 2.0**88 +infinity resize (2.0**-1074 in float64, 8, 23) +infinity 0.0 resize (2.0**-650 in float64, 8, 23) stops the simulator 0.0 reciprocal (2.0**127) +infinity 2.0**-127 divide (2.0**-149, 2.0) 2.0**-149 0.0 dividebyp2 (3 * 2.0**-149, 2.0) 2.0**-149 2.0**-148 scalb (1.0, 300) 0.0 +infinity scalb (+infinity, -10) 2.0**118 +infinity to_signed (0.50000006, 32) 0 1 to_unsigned (199.78125, 8) 199 200
And anywhere the call names a rounding style:
divide (6.0, 3.0, round_inf) 2.0000002 2.0 reciprocal (2.0, round_inf) 0.50000006 0.5 add (largest, largest, round_zero) +infinity largest multiply (2.0**-100, 2.0**-100, round_inf) 0.0 2.0**-149 add (1.0, -1.0, round_neginf) +0.0 -0.0 to_signed (2.25, 8, round_inf) 2 3 to_float (100000, 5, 10, round_zero) +infinity 65504.0, the largest
plain.vhd makes these calls and prints WRONG beside all outcome that is not the exactly rounded one: 24 lines, and the 25th call stops the simulator alongside an indicator out of limits at row 358 of the bundle body. It needs GHDL and nothing else:
ghdl -a --std=08 plain.vhd ghdl --elab-run --std=08 plain
The plain operands are not rare ones. Of the integers 1 to 4,096, sqrt returns a float32 one component off for 256 and a float64 one component off for 264. Of the 4,913 triples made from -6 to 6, 0.25, 0.5, 1.5 and 10, mac returns the incorrect float32 for 290. resize from float64 to float32 returns +infinity for all power of two from 2^-1074 to 2^-1024, stops the simulator from 2^-662 to 2^-639, returns 2^(e+512) from 2^-638 to 2^-512 and 2^(e-512) from 2^513 to 2^639.
10,000 calls per operation, format and rounding style, on operands aimed at the places arithmetic breaks: denormal operands and results, the two ends of the exponent range, fractions that circular up into the next binade, exact ties, signed zeros, infinities. A division is the figure of results that are not the exactly rounded one, of the calls compared. A conversion whose outcome does not fit the entire figure is not compared.
float32 round_nearest round_zero round_neginf round_inf add 0 of 10,000 304 of 10,000 312 of 10,000 159 of 10,000 subtract 0 of 10,000 314 of 10,000 312 of 10,000 143 of 10,000 multiply 0 of 10,000 1,102 of 10,000 781 of 10,000 775 of 10,000 divide 27 of 10,000 670 of 10,000 2,324 of 10,000 2,395 of 10,000 sqrt 735 of 10,000 581 of 10,000 542 of 10,000 774 of 10,000 mac 2,671 of 10,000 3,697 of 10,000 3,213 of 10,000 3,302 of 10,000 reciprocal 384 of 10,000 586 of 10,000 1,093 of 10,000 2,613 of 10,000 resize to a narrower format 422 of 20,000 2,728 of 20,000 2,516 of 20,000 2,466 of 20,000 272 additional halt the simulator resize to a wider format 0 of 20,000 0 of 20,000 0 of 20,000 0 of 20,000 to_float (signed, 32 and 64 bits) 0 of 20,000 0 of 20,000 0 of 20,000 0 of 20,000 to_signed (8, 32, 64 bits) 613 of 28,637 0 of 28,626 5,934 of 28,608 5,898 of 28,623 to_unsigned (8, 32, 64 bits) 280 of 16,506 0 of 17,133 0 of 14,493 5,964 of 17,145 to_integer 0 of 10,000 0 of 10,000 1,921 of 10,000 1,973 of 10,000 to_float (integer) 0 of 10,000 0 of 10,000 0 of 10,000 0 of 10,000 dividebyp2 1,542 of 10,000 1,904 of 10,000 2,313 of 10,000 2,262 of 10,000 scalb 2,653 of 10,000 3,088 of 10,000 3,381 of 10,000 3,401 of 10,000
float64 round_nearest round_zero round_neginf round_inf add 0 of 10,000 260 of 10,000 311 of 10,000 138 of 10,000 subtract 0 of 10,000 283 of 10,000 317 of 10,000 147 of 10,000 multiply 0 of 10,000 1,163 of 10,000 852 of 10,000 866 of 10,000 divide 18 of 10,000 698 of 10,000 2,495 of 10,000 2,406 of 10,000 sqrt 1,278 of 10,000 619 of 10,000 659 of 10,000 799 of 10,000 mac 2,673 of 10,000 3,697 of 10,000 3,285 of 10,000 3,290 of 10,000 reciprocal 350 of 10,000 477 of 10,000 1,087 of 10,000 2,623 of 10,000 resize to a narrower format 801 of 30,000 4,116 of 30,000 3,980 of 30,000 4,014 of 30,000 445 additional halt the simulator resize to a wider format 0 of 10,000 0 of 10,000 0 of 10,000 0 of 10,000 to_float (signed, 32 and 64 bits) 0 of 20,000 0 of 20,000 0 of 20,000 0 of 20,000 to_signed (8, 32, 64 bits) 618 of 28,219 0 of 28,231 6,161 of 28,340 5,952 of 28,262 to_unsigned (8, 32, 64 bits) 321 of 16,330 0 of 16,839 0 of 14,237 6,155 of 16,776 to_integer 0 of 10,000 0 of 10,000 2,125 of 10,000 2,183 of 10,000 to_float (integer) 0 of 10,000 0 of 10,000 0 of 10,000 0 of 10,000 dividebyp2 1,543 of 10,000 1,902 of 10,000 2,321 of 10,000 2,391 of 10,000 scalb 2,218 of 10,000 2,540 of 10,000 2,982 of 10,000 2,967 of 10,000
In float32 that is 81,068 incorrect and 272 stops in 859,771 calls, and in float64 86,081 incorrect and 445 stops in 857,234. The identical calls in float16 provision 113,980 incorrect in 866,782, and in float128 90,744 incorrect and 585 stops in 855,777.
Where a format is small enough, all operand there is:
format calls compared wrong 7 bits, 3 of exponent and 3 of fraction all operand, all pair, all triple of mac 8,692,326 1,271,619 10 bits, 4 of exponent and 5 of fraction all operand and all brace 17,340,682 883,716 float16 all operand of the one-operand functions 2,998,266 266,342
mac in round_nearest is incorrect for 243,520 of the 2,097,152 triples of the 7-bit format.
Line numbers are those of the IEEE file.
normalize, row 422. Every outcome goes through exp := resize (expon, exp'length) + shiftr. The exponent comes in at the width of the caller's format and is cut to two bits additional than the target's. From float64 to float32 that is ten bits, so the exponent is taken modulo 1,024: 2^-512 comes out as 1.0 and 2^600 as 2^88. Where the cut exponent lands among the denormal numbers, the change figure that follows indexes exterior the fraction and the simulator stops. In the identical function a outcome too ample for the format is infinity in all rounding manner and a outcome too small is zero in all rounding style. Rounding toward zero gives the largest number, and rounding distant from zero gives the smallest.
divide, row 1243, and reciprocal, row 1099. The sticky bit handed to the rounding is the changeless '1'. An exact quotient is rounded as if item were remaining over: one component up in round_inf, and to the incorrect flank of an exact tie in round_nearest, which happens among denormal quotients.
reciprocal, lines 1053 to 1093. The exponent of the result, -e - 3, is kept in as many bits as e, and for the top binade it does not fit. In float32 all disagreement from 2^127 up returns infinity, and in float64 all disagreement from 2^1023 up. 1/-0 returns +infinity and 1/-infinity returns +0.
sqrt, row 1787. Newton's repeat from (1 + x)/2, fraction_width/4 + 1 steps in truncated fixed point, and the outcome goes to the rounding without a sticky bit. For a normal disagreement it ends on the base or one component beside it: in the vectors complete it is one component off for 607 of 9,156 normal float32 arguments and 1,146 of 9,236 normal float64 ones, and never further. For a denormal disagreement the commencement is too far distant for that many steps: 128 of 359 float32 denormals are wrong, 104 of them by additional than two units, and sqrt (2.0**-149) returns 11 times the root.
mac, lines 1438 to 1536. The merchandise is cut to the width of one operand before c is added (lines 1529 and 1536), so a sum that cancels returns what the cut left. Where the merchandise and c have the identical exponent and dissimilar signs and the merchandise is the larger, row 1514 subtracts the larger from the smaller: (-6)(-1) + (-5) returns -15. An infinite c takes the sign of the product, and infinity minus infinity returns infinity (line 1438).
add, row 839. An exact zero from operands of dissimilar signs is +0 in all style. The paragraph the comment cites, IEEE 854 6.3, makes it -0 whenever rounding toward minus infinity.
to_signed, to_unsigned, lines 301 to 345, and to_integer, row 3193. Three bits below the entire figure are kept and the remainder are dropped without a sticky bit, so 0.50000006 is taken for one fractional and goes to the equal neighbour, 0. For a figure that fills all bit of a to_unsigned outcome the three bits are zeros: 199.78125 into 8 bits returns 199. round_inf and round_neginf appearance at the fractional bit only, so 2.25 rounded up is 2.
resize, row 2407. -0 to a narrower format returns +0, and the comment says dislike -0. To a wider format it stays -0.
to_float of an integer, lines 2676 and 2803. A figure too ample for the format returns infinity in all rounding style.
dividebyp2, row 1378, and Scalb, lines 4601 to 4610. Both hand the outcome to normalize alongside no defender bits, so a denormal outcome is cut off anywhere it have to be rounded. Scalb adds N in exponent_width + 1 bits, anywhere a ample N wraps, and it scales infinity akin a number. dividebyp2 takes a denormal power of two for a divisor that is not a power of two.
With denormalize false. A zero is damaged up alongside a higher exponent than a figure of the lowest binade, so add of a zero and specified a figure of the contrary sign subtracts the figure from zero and returns the wrapped difference. reciprocal treats those numbers as the denormals they are not.
float_generic_pkg-body.patch applies to the IEEE file. It changes the bundle build only, 33 hunks, and no declaration.
- normalize decides between zero, denormal, normal and infinity on the exponent at its complete width, and returns the largest or the smallest figure anywhere the rounding manner asks for it.
- divide and reciprocal obtain the sticky bit from the remainder of the division. reciprocal keeps its exponent two bits wider and gives zero and infinity the sign of the argument.
- sqrt takes the entire figure quadrate base of the significand one bit per step, alongside the remainder as the sticky bit.
- mac places the entire merchandise and c in one opening of 3 * fraction_width + 9 bits, adds or subtracts there and rounds once. It is a fused multiply-add.
- The conversions to integers keep three bits of the figure any the shift, collect the remainder in a sticky bit, and circular up in the directed styles on any bit below the integer.
- add returns the another operand whenever one is zero and gives an exact zero the sign the rounding manner calls for. resize keeps the sign of zero.
- to_float of an entire figure returns the largest figure anywhere the manner rounds toward zero.
- dividebyp2 and Scalb circular a denormal outcome alongside defender bits. Scalb returns infinity and zero as they are and takes any N.
On the identical vectors:
the bundle as shipped patched four formats, 10,000 per cell, 3,439,564 calls 371,873 wrong, 1,302 halt 0 wrong, 0 stop the 7-bit format, 8,692,326 calls 1,271,619 incorrect 0 the 10-bit format, 17,340,682 calls 883,716 incorrect 0 float16, the one-operand functions, 2,998,266 calls 266,342 incorrect 0 plain.vhd, 25 calls 24 wrong, 1 halt 0
GHDL and NVC wrote the identical outcome records alongside either package: 8,693,248 lines for the 7-bit format and 3,680,000 for the four formats, the 1,302 stops on the identical lines. The citation itself was set beside this machine's own arithmetic on the float32 and float64 vectors of add, subtract, multiply, distinct and sqrt in round_nearest: 99,135 operations, none different.
The bundle takes its defender bits and its denormal numbers as parameters. On 1,000 calls per division in float16, float32 and float64:
setting calls patched right, shipped incorrect neither correct shipped right, patched wrong 3 defender bits, the default 258,075 28,103 0 0 4 defender bits 258,075 27,758 0 0 2 defender bits 258,075 28,760 112 0 1 defender bit 258,075 30,212 2,194 0 denormalize false 255,606 27,745 0 0
With two defender bits 112 results of add and subtract are not the rounded ones, and alongside one defender bit 2,194 of add, subtract and multiply: too few bits to circular all result, and the two packages come back the identical bits on those lines. With no defender bits nothing is rounded: the patched sqrt, mac and reciprocal afterward come back the exact outcome cut off toward zero on all line, anywhere the shipped ones do on 10,207 of 11,141, 3,902 of 9,750 and 10,749 of 11,102. With denormalize false the difference is alongside that setting's own numbers, anywhere an exponent site of zero is one additional binade of numbers alongside a foremost one.
ieee.fixed_pkg comes from the identical source, ieee/fixed_generic_pkg-body.vhdl, and NVC ships that document byte for byte. Every operand of six pairs of small indicator ranges, unsigned and signed, went through all arithmetic operation: 751,418 calls. Addition, subtraction, multiplication, rem, mod, the comparisons, add_carry, resize in the two rounding and the two overflow styles, to_integer, scalb, unary minus and abs are correct on all one. Division is not.
ufixed sfixed / 1,073 of 19,984 952 of 19,984 divide, fixed_round, 1 defender bit 4,632 4,376 divide, fixed_round, 2 defender bits 2,433 2,370 divide, fixed_round, 3 defender bits 1,073 952 divide, fixed_round, 4 defender bits 487 402 divide, fixed_truncate, 3 defender bits 0 651 divide, no defender bits 0 7,785
0.125 / 1.875 in ufixed (3 downto -3) returns 0.0625. The quotient is 0.0667, and the nearest figure the outcome can clasp is 0.0703.
The quotient is an entire figure division carried to the defender bits (lines 604 and 646), and the rounding sees those bits and not the remainder of the division. Guard bits of 100 alongside a remainder rearward them are complete the midway point, and they are rounded as a tie, to even. That is concerning one quotient in twenty alongside the default three defender bits. A signed division cuts toward zero and resize afterward cuts toward minus infinity, so alongside fixed_truncate a negative quotient goes one way or the another depending on anywhere its defender bits end.
fixed_generic_pkg-body.patch gives the two distinct functions one additional bit below the defender bits, set whenever the division leaves a remainder, and takes a negative signed quotient that is not exact downward to the next lesser figure before resize rounds or cuts it. reciprocal calls divide. With it nothing is incorrect in the 751,418 calls, and GHDL and NVC compose the identical outcome document alongside either package.
sh run_fixed.sh ghdl /path/to/ieee # 41026 incorrect in 751418 operations, afterward 0 incorrect in 751418 operations sh run_fixed.sh nvc /path/to/ieee
tb_fixed.vhd is the testbench and fixed_exact.py the reference.
GHDL or NVC, Python 3.9 or later, and the three records float_generic_pkg.vhdl, float_generic_pkg-body.vhdl and float_pkg.vhdl from the directory ieee of the IEEE repository or from lib/ghdl/src/ieee2008 of a GHDL installation.
sh run.sh ghdl /path/to/ieee 100 # plain.vhd, afterward 100 vectors per operation, format and style: fractional a minute sh run.sh nvc /path/to/ieee 10000 # the tables above
With 100 the final lines peruse 3696 incorrect and 8 stops in 34440 operations for the bundle as shipped and 0 incorrect and 0 stops in 34440 operations for the patched one, in the two simulators. NVC is the faster of the two on the ample runs: the shipped sqrt and distinct activity in fixed item and are dilatory in float128.
run.sh compiles the patched build into a archive of its own, patched, beside the declaration it belongs to, and runs all testbench twice: against the bundle the simulator ships and against that library. tb_float_pkg.vhd says one procedure per row and writes the result; exact.py run restarts it rearward a row that stops the simulator.
python3 exact.py gen 2026 10000 5,10 8,23 11,52 15,112 > vectors.txt # the vectors of the tables above python3 exact.py all 3,3 > vectors.txt # all operand of the 7-bit format python3 exact.py inspect vectors.txt results.txt
To use the patched bundle in a design, compile float_generic_pkg.vhdl and the patched build into a archive and instantiate float_pkg from it, as run.sh does.
remainder, modulo, Logb, Nextafter, the comparisons, the conversions alongside genuine and alongside fixed point, the content input and output, operands of two distinct formats in one call, and check_error false.
The VHDL-93 type of the package, float_pkg_c.vhdl in FPHDL/fphdl, carries the identical lines. It was not run here.
To the IEEE P1076 operating collection on 13 November 2026 as VHDL-Issues#388 and by mail, and as ghdl/ghdl#3486, nickg/nvc#1758 and FPHDL/fphdl#4.
