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

[Craig Bruce]

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".

If you want, you could atof() the re-generated string back to a double
and check to see if it is exactly equal to the given value.  If it is,
then the 15-digit value is definitely the correct value.  If not, then
it means that the .0000 and .0001 representations are distinct and the
"ugly" .0001 representation can be used instead, though as noted below,
this presumes that the source value was correct to begin with.  This step
is optional.

This approach gives you nearly full double precision all of the time
without the eyesores.

Here is a C function that implements what I am talking about (plus it
takes precision as a parameter).  (It uses a 'bool' type with values
'true' and 'false'.  'cw_traceN' are debugging macros.)  It is licensed
under the LGPL.

/*------------------------------------------------------------------------*\
 *  NAME:
 *     CwStr_GenDoublePrec() - generate a 'double' to a string buffer
 *  DESCRIPTION:
 *     Generates a double-precision floating-point number to a decimal
 *     text representation similar to printf's "%.g" into a string buffer
 *     with the given number of significant digits, which is clipped
 *     to be between 1 and 20.  If 16 significant digits is selected,
 *     then error control is applied to favor rounder decimal numbers
 *     when the 15th and 16th digits have the pattern '01' or '99'.
 *     Since 'double's only have around 16 significant decimal digits
 *     of precision, selecting more will only show artifacts in the IEEE
 *     'double' encoding, which may be desirable in some circumstances.
 *  NOTE:
 *     It would be desirable to verify that the corrected value has
 *     exactly the same 'double' bit pattern as the given value and to
 *     undo the correction if this is not true.  This would mean that
 *     the given value is numerically distinct from the rounded value.
 *     Unfortunately, there are faulty sources of double-precision
 *     values, such as the Oracle database.  E.g., given the Oracle
 *     decimal NUMBER value 4.877, it generates the wrong 'double' value,
 *     4.8770000000000007 instead of the correct 4.8769999999999998.
 *     Verification is not practical because of faulty sources like this.
 *  ARGUMENTS:
 *     strBuf - string buffer to generate number into, at least [nDigits+8]
 *     number - number to print
 *     nDigits - number of significant digits to generate
 *  RETURNS:
 *     strBuf - pointer to given string buffer
 *  ERRORS:
 *     (no errors are possible)
\*------------------------------------------------------------------------*/

char *CwStr_GenDoublePrec( char *strBuf, double number, long nDigits )
{
    static const char *printfFmts[] = { NULL, "%.1g", "%.2g", "%.3g", "%.4g",
        "%.5g", "%.6g", "%.7g", "%.8g", "%.9g", "%.10g", "%.11g", "%.12g",
        "%.13g", "%.14g", "%.15g", "%.16g", "%.17g", "%.18g", "%.19g", "%.20g"};
    const char *printfFmt;
    bool amAfterDecimal, isZero;
    char c, numChar15;
    long nScanDigits, pos;
#if 0
    double verifyNum;
#endif

    cw_trace2(("CwStr_GenDoublePrec(number=%.16g, nDigits="COUNT_FMT")\n",
                number, nDigits ));

    /** generate raw string **/
    if (nDigits < 1) nDigits = 1;
    if (nDigits > 20) nDigits = 20;
    printfFmt = printfFmts[nDigits];
    sprintf( strBuf, printfFmt, number );

    /** apply error control **/
    if (nDigits == 16) {
        /** scan number for pattern of digits 15 & 16 **/
        cw_trace3(("checking if number needs to be corrected\n"));
        amAfterDecimal = false;
        nScanDigits = 0;
        c = 'x';
        numChar15 = 'x';
        isZero = true;
        for (pos=0; strBuf[pos] != '\0'; pos++) {
            c = strBuf[pos];
            if (c >= '0' && c <= '9') {
                isZero = (isZero && c == '0');
                if (!isZero) {
                    nScanDigits += 1;
                    if (nScanDigits == 15) {
                        numChar15 = c;
                    } else if (nScanDigits == 16) {
                        break;
                    }
                }
            } else if (c == '.') {
                amAfterDecimal = true;
            } else if (c == 'e') {
                break;
            }
        }
        cw_trace3(("scan: nScanDigits=%d, amAfterDecimal=%s, numChar15='%c', "
            "c='%c'\n", nScanDigits, CwStr_DecodeBool(amAfterDecimal),
             numChar15, c ));

        /** correct number if necessary **/
        if (nScanDigits == 16 && amAfterDecimal
                && ((numChar15 == '0' && c == '1')
                    || (numChar15 == '9' && c == '9'))) {
            /** correct '01' or '99' pattern **/
            cw_trace3(("correcting '...%c%c' pattern\n", numChar15, c ));
            sprintf( strBuf, "%.15g", number );

#if 0
            /** verify correction **/
            cw_trace3(("verifying corrected number \"%s\"\n", strBuf ));
            verifyNum = atof(strBuf);
            if (verifyNum == number) {  /* exact match */
                cw_trace3(("verification confirmed; number corrected: "
                    "num=%.20g, vfy=%.20g, strBuf=\"%s\"\n", number,
                     verifyNum, strBuf ));
            } else {
                sprintf( strBuf, "%.16g", number );
                cw_trace3(("verification failed; number restored: "
                    "num=%.20g, vfy=%.20g, strBuf=\"%s\"\n", number,
                     verifyNum, strBuf ));
            }
#endif
        }
    }
    return( strBuf );
}

--------------------------+----------------------+--------------------------
Dr. Craig S. Bruce        | Ph 819-771-8303 x205 |             CubeWerx Inc.
Senior Software Developer |   Fax 819-771-8388   |  Gatineau, Québec, Canada
csbruce at cubewerx.com      |  http://csbruce.com/ |  http://www.cubewerx.com/
--------------------------+----------------------+--------------------------
                "I am both amused and annoyed that you think
            I should be less stubborn than you are." -- House, MD