[metapost] Units unaccuracy with numbersystem set to double
Nelson H. F. Beebe
beebe at math.utah.edu
Thu Apr 3 16:45:47 CEST 2014
Franck Pastor <franck.pastor at mac.com> writes today in
response to my earlier posting
>> ...
>> > A more accurate conversion from PostScript big points to cm is this:
>> >
>> > hocd128> 72 / 2.54
>> > 28.346_456_692_913_385_826_771_653_543_307_09
>> >
>> > Now that Metapost supplies arithmetic of higher precision, perhaps
>> > input files and source code should be examined for truncated
>> > floating-point constants, and those constants replaced by expressions,
>> > such as the exact 72 / 2.54 above, that could then be evaluated at run
>> > time in the prevailing precision, which might be higher than that in
>> > classic METAFONT.
>> >
>>
>> 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.
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
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:
@InCollection{Knuth:1990:SPW,
author = "Donald E. Knuth",
title = "A Simple Program Whose Proof Isn't",
crossref = "Feijen:1990:BOB",
chapter = "27",
pages = "233--242",
year = "1990",
bibdate = "Mon Feb 03 07:07:55 2003",
bibsource = "http://www.math.utah.edu/pub/tex/bib/fparith.bib",
note = "This paper discusses the algorithm used in {\TeX} for
converting between decimal and scaled fixed-point
binary values, and for guaranteeing a minimum number of
digits in the decimal representation. See also
\cite{Clinger:1990:HRF,Clinger:2004:RHR} for decimal to
binary conversion,
\cite{Steele:1990:HPF,Steele:2004:RHP} for binary to
decimal conversion, and \cite{Gries:1990:BDO} for an
alternate proof of Knuth's algorithm.",
acknowledgement = ack-nhfb,
keywords = "decimal floating-point arithmetic",
}
@String{pub-SV = "Spring{\-}er-Ver{\-}lag"}
@String{pub-SV:adr = "Berlin, Germany~/ Heidelberg,
Germany~/ London, UK~/ etc."}
@Book{Feijen:1990:BOB,
editor = "W. H. J. Feijen and A. J. M. {van Gasteren} and D.
Gries and J. Misra",
booktitle = "Beauty is our business: a birthday salute to {Edsger
W. Dijkstra}",
title = "Beauty is our business: a birthday salute to {Edsger
W. Dijkstra}",
publisher = pub-SV,
address = pub-SV:adr,
pages = "xix + 453",
year = "1990",
ISBN = "0-387-97299-4",
ISBN-13 = "978-0-387-97299-2",
LCCN = "QA76 .B326 1990",
bibdate = "Thu Mar 24 09:27:40 1994",
bibsource = "http://www.math.utah.edu/pub/tex/bib/fparith.bib",
acknowledgement = ack-nhfb,
}
The references in the \cite{} macros are all in fparith.bib at the
indicated URL (change .bib to .html for a live hyperlink display in a
Web browser).
-------------------------------------------------------------------------------
- Nelson H. F. Beebe Tel: +1 801 581 5254 -
- University of Utah FAX: +1 801 581 4148 -
- Department of Mathematics, 110 LCB Internet e-mail: beebe at math.utah.edu -
- 155 S 1400 E RM 233 beebe at acm.org beebe at computer.org -
- Salt Lake City, UT 84112-0090, USA URL: http://www.math.utah.edu/~beebe/ -
-------------------------------------------------------------------------------
More information about the metapost
mailing list