[Proj] Numerical Precision in Proj, WKT, and everywhere

[Karney, Charles]

What you want is a "perfect floating-point print" routine.  This takes a
binary floating point number and prints a decimal representation with the
following properties:

* It's a correctly rounded version of the exact decimal number.

* When read by a perfect floating-point read routine, the original
  binary representation is recovered.

* It's the shortest representation satisfying these properties.

This problem was tackled by White and Steele in the implementation of
the MacLisp PRINT function in the 1970s.  The algorithms are given in:

    G. L. Steele and J. L. White
    How to print floating-point numbers accurately
    ACM SIGPLAN 1990
    http://dx.doi.org/10.1145/93542.93559
    See also:
    http://dx.doi.org/10.1145/989393.989431

It turns out that the "perfect read" routine is not trivial, see

    W. D. Clinger
    How to read floating point numbers accurately
    ACM SIGPLAN 1990
    http://dx.doi.org/10.1145/93542.93557

--
Charles Karney <ckarney at sarnoff.com>
Sarnoff Corporation, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662
________________________________________
From: proj-bounces at lists.maptools.org [proj-bounces at lists.maptools.org] On Behalf Of Craig Bruce [csbruce at cubewerx.com]
Sent: Wednesday, March 31, 2010 19:36
To: proj at lists.maptools.org
Subject: [Proj] Numerical Precision in Proj, WKT, and everywhere

I was reading the following web page and its discussion of numerical
precision written by Frank Warmerdam:

http://home.gdal.org/projects/opengis/wktproblems.html

I agree with what he said, but printf's "%.16g" does have a problem which
manifests itself frequently, such as in the "epsg" CRS definitions in
proj-4.6.0, e.g.:

# Anguilla 1957 / British West Indies Grid
<2000> +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +units=m +no_defs  <>

The "+k" parameter value is an eyesore.  I think a better approach is to
generate the "%.16g" number as indicated, but then scan its characters for
the ugly pattern of '01' or '99' in the 15th and 16th significant digits.
If this pattern is found, then the number should be generated again using
"%.15g".

**SNIP**