[metapost] Units unaccuracy with numbersystem set to double

[Dan Luecking]

On Thu Apr 3, 2014, at 16:45 CEST Nelson H. F. Beebe <beebe at math.utah.edu>
wrote:


>The result of computing 72 / 2.54 is slightly less accurate: it is
>about 5 ulps in error:
>
>         hocd128> ((28.34653 - ecm)/ecm) * 2**21
>                  5.423_468_088_888_888_888_888_888_888_638_394
>
>That happens because of the conversion of the denominator to
>a fixed-point value:
>
>     hocd128> floor(2.54 * 2**16)
>              166_461
>
>     hocd128> 166_461 / 2**16
>              2.539_993_286_132_812_5
>
>     hocd128> 72 / 2.539_993_286_132_812_5
>              28.346_531_620_019_103_573_810_081_640_744_68
>
>That is why 72 / 2.54 evaluates to 28.34653.

Some minor corrections: Knuth claims that decimal numbers
are converted to the closest multiple of $2^{-16}$. Thus
"floor(2.54 * 2**16)" should rather be "round(2.54 * 2**16)".
Since the decimal part is .44 in this case, the two formulas
produce the same result.

Also, the last two computations in the scaled system are
actually computed as (essentially) 72*(2**16)/166_461 and
this produces 1_857_718/(2**16). This last fraction is
(exactly) 28.346_527_099_609_375, but will be converted to
decimal as the value 28.34653 as that is the shortest decimal
that reproduces the internal 1857718/(2**16).


Dan


Daniel H. Luecking
Department of Mathematical Sciences
Fayetteville, Arkansas
http://www-cs-faculty.stanford.edu/~knuth/iaq.html