[metapost] Units unaccuracy with numbersystem set to double

[Franck Pastor]

Le 3 avr. 2014 à 16:45, Nelson H. F. Beebe a écrit :

> Franck Pastor <franck.pastor at mac.com> writes today in
>>> 
>>> It seems it is s unfortunately more complicated than that. With the
>>> default arithmetics of MetaPost, 72/2.54 gives a result noticeably
>>> less accurate than the constant value of "cm" to be found in
>>> plain.mp. For comparison:
>>> 
>>> *show 72/2.54;
>>>>> 28.34653
>>> *show cm;
>>>>> 28.34645
>>> ...
> 
> That happens because METAFONT and MetaPost by default work with scaled
> binary arithmetic, where a 32-bit word is divided into a sign bit, an
> overflow bit, a 14-bit integer part, and a 16-bit fractional part.
> 
> 	hocd128> ecm = 72 / 2.54
> 
> 	hocd128> ecm * 2**16
> 	         1_857_713.385_826_771_653_543_307_086_614_173
> 
> 	hocd128> floor(ecm * 2**16) / 2**16
> 	         28.346_450_805_664_062_5
> 
> It is the limited precision of the fraction that is responsible for
> the printed value of cm.
> 
> Because 5 bits are needed for the integer part of cm, we effectively
> have a 21-bit approximation to the correct value.  The two decimal
> values of ecm and cm are closer than half a unit in the last place:
> 
> 	hocd128> ((28.34645 - ecm) / ecm) * 2**21
> 	         -0.495_160_888_888_888_888_888_888_889_139_381_9
> 
> Thus, either one is acceptable for the default fixed-point arithmetic.
> We can demonstrate that in Metapost:
> 
> 	*xcm := 28.346456692913;
> 
> 	*show xcm;
> 	>> 28.34645
> 
> 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.


Impressive! And very clear. I couldn't have understood that by myself! 
> 
> The point of my posting is that such things confuse users, and now
> that Metapost offers alternative higher-precision arithmetics, it
> would be best to avoid use of inexact fractional constants in Metapost
> programs.  Thus, 72 / 2.54 (where 2.54 is EXACT by an international
> standard) might be the preferred exact way of expressing the value of
> cm (which is really big_points_per_cm), because it is an exact
> rational number that can be then converted to whatever the working
> precision is in the current run of the program.
> 
> I would instead write cm as 7200 / 254, so that no unnecessary
> rounding errors are introduced.  However, the value 7200 then exceeds
> the 14-bit limit on the integer part.  By eliminating a common factor
> of 2 in the ratio, we get 3600 / 127, a value that IS acceptable to
> Metapost:
> 
> 	*x := 3600 / 127;
> 
> 	*show x;
> 	>> 28.34645


It looks like it is the way to go for the "cm" definition. And "mm" would be ideally defined as 360/127 or cm/10 if I understand correctly…

> 
> I dug into METAFONT's plain.mf and modes.mf, and found that cm is set
> in a round-about way:
> 
> 	cm := pixels_per_inch / 2.54;
> 
> where pixels_per_inch comes from a device-specific setting.  However,
> in Metapost, mfplain.mp simply has the assignments
> 
> 	cm #= 28.34645
> 	in #= 72
> 
> that correspond to units of big points.
> 
> Don Knuth wrote a lovely paper with a memorable title on the
> fixed-point binary arithmetic of TeX and Metafont, and the problems of
> representing numbers in decimal such that exact round-trip conversion
> (decimal -> binary -> decimal) is guaranteed:

It's probably way above my head :-) I think I will first re-read my old books about numerical analysis and most of all their chapters about rounding errors. And then I'll go back to this article. I'm more and more interested in this subject.

Many thanks,

Franck Pastor