Support for high precision arithmetic

Back to Auxiliary latitudes. Forward to Change log. Up to Contents.

One of the goals with the algorithms in GeographicLib is to deliver accuracy close to the limits for double precision. In order to develop such algorithms it is very useful to be have accurate test data. For this purpose, I used Maxima's bfloat capability, which support arbitrary precision floating point arithmetic. As of version 1.37, such high-precision test data can be generated directly by GeographicLib by compiling it with GEOGRAPHICLIB_PRECISION equal to 4 or 5.

Here's what you should know:

• This is mainly for use for algorithm developers. It's not recommended for installation for all users on a system.
• Configuring with -D GEOGRAPHICLIB_PRECISION=4 gives quad precision (113-bit precision) via boost::multiprecision::float128; this requires:
  • Boost, version 1.64 or later,
  • the quadmath library (the package names are libquadmath and libquadmath-devel),
  • the use of g++.
• Configuring with -D GEOGRAPHICLIB_PRECISION=5 gives arbitrary precision via mpfr::mpreal; this requires:
  • MPFR, version 3.0 or later,
  • MPFR C++ (version 3.6.9, dated 2022-01-18, or later; version 3.6.9 also requires the fixes given in pull requests #15),
  • a compiler which supports the explicit cast operator (e.g., g++ 4.5 or later, Visual Studio 12 2013 or later).
• MPFR, MPFR C++, and Boost all come with their own licenses. Be sure to respect these.
• The indicated precision is used for all floating point arithmetic. Thus, you can't compare the results of different precisions within a single invocation of a program. Instead, you can create a file of accurate test data at a high precision and use this to test the algorithms at double precision.
• With MPFR, the precision should be set (using Utility::set_digits) just once before any other GeographicLib routines are called. Calling this function after other GeographicLib routines will lead to inconsistent results (because the precision of some constants like Math::pi() is set when the functions are first called).
• All the Utility programs call Utility::set_digits() (with no arguments). This causes the precision (in bits) to be determined by the GEOGRAPHICLIB_DIGITS environment variable. If this is not defined the precision is set to 256 bits (about 77 decimal digits).
• The accuracy of most calculations should increase as the precision increases (and typically only a few bits of accuracy should be lost). We can distinguish 4 sources of error:
  • Round-off errors; these are reliably reduced when the precision is increased. For the most part, the algorithms used by GeographicLib are carefully tuned to minimize round-off errors, so that only a few bits of accuracy are lost.
  • Convergence errors due to not iterating certain algorithms to convergence. However, all iterative methods used by GeographicLib converge quadratically (the number of correct digits doubles on each iteration) so that full convergence is obtained for "reasonable" precisions (no more than, say, 100 decimal digits or about 340 bits). An exception is thrown if the convergence criterion is not met when using high precision arithmetic.
  • Truncation errors. Some classes (namely, Geodesic and TransverseMercator) use series expansion to approximate the true solution. Additional terms in the series are used for high precision, however there's always a finite truncation error which varies as some power of the flattening. On the other hand, GeodesicExact and TransverseMercatorExact are somewhat slower classes offering the same functionality implemented with EllipticFunction. These classes provide arbitrary accuracy. (However, a caveat is that the evaluation of the area in GeodesicExact still uses a series (albeit of considerably higher order). So the area calculations are always have a finite truncation error.)
  • Quantization errors. Geoid, GravityModel, and MagneticModel all depend on external data files. The coefficient files for GravityModel and MagneticModel store the coefficients as IEEE doubles (and perhaps these coefficients can be regarded as exact). However, with Geoid, the data files for the geoid heights are quantized at 3mm leading to an irreducible ±1.5mm quantization error. On the other hand, all the physical constants used by GeographicLib, e.g., the flattening of the WGS84 ellipsoid, are evaluated as exact decimal numbers.
• Where might high accuracy be important?
  • checking the truncation error of series approximations;
  • checking for excessive round-off errors (typically due to subtraction);
  • checking the round-off error in computing areas of many-sided polygons;
  • checking the summation of high order spherical harmonic expansions (where underflow and overflow may also be a problem).
• Because only a tiny number of people will be interested in using this facility:
  • the cmake support for the required libraries is rudimentary;
  • however geographiclib-config.cmake does export GEOGRAPHICLIB_PRECISION and GeographicLib_HIGHPREC_LIBRARIES, the libraries providing the support for high-precision arithmetic;
  • support for the C++11 mathematical functions and the explicit cast operator is required;
  • quad precision is only available on Linux;
  • mpfr has been mostly tested on Linux (but it works on Windows with Visual Studio 12 and MacOS too).

The following steps needed to be taken

• Phase 1, make sure you can switch easily between double, float, and long double.
  • use #include <cmath> instead of #include <math.h>;
  • use, e.g., std::sqrt instead of sqrt in header files (similarly for sin, cos, atan2, etc.);
  • use using namespace std; and plain sqrt, etc., in code files;
  • express all convergence criteria in terms of

    numeric_limits<double>::epsilon() 

    etc., instead of using "magic constants", such as 1.0e-15;
  • use typedef double real; and replace all occurrences of double by real;
  • write all literals by, e.g., real(0.5). Some constants might need the L suffix, e.g., real f = 1/real(298.257223563L) (but see below);
  • Change the typedef of real to float or long double, compile, and test. In this way, the library can be run with any of the three basic floating point types.
  • If you want to run the library with multiple floating point types within a single executable, then make all your classes take a template parameter specifying the floating-point type and instantiate your classes with the floating-point types that you plan to use. I did not take this approach with GeographicLib because double precision is suitable for the vast majority of applications and turning all the classes into templates classes would end up needlessly complicating (and bloating) the library.
• Phase 2, changes to support arbitrary, but fixed, precision
  • Use, e.g.,

    typedef boost::multiprecision::float128 real; 

  • Change std::sqrt(...), etc. to

    using std::sqrt;
    sqrt(...) 

    (but note that std::max can stay). It's only necessary to do this in header files (code files already have using namespace std;).
  • In the case of boost's multiprecision numbers, the C++11 mathematical functions need special treatment, see Math.hpp.
  • If necessary, use series with additional terms to improve the accuracy.
  • Replace slowly converging root finding methods with rapidly converging methods. In particular, the simple iterative method to determine the flattening from the dynamical form factor in NormalGravity converged too slowly; this was replaced by Newton's method.
  • If necessary, increase the maximum allowed iteration count in root finding loops. Also throw an exception of the maximum iteration count is exceeded.
  • Write literal constants in a way that works for any precision, e.g.,

    real f = 1/( real(298257223563LL) / 1000000000 ); 

    [Note that

    real f = 1/( 298 + real(257223563) / 1000000000 ); 

    and 1/real(298.257223563L) are susceptible to double rounding errors. We normally want to avoid such errors when real is a double.]
  • For arbitrary constants, you might have to resort to macros

    #if   GEOGRAPHICLIB_PRECISION == 1
    #define REAL(x) x##F
    #elif GEOGRAPHICLIB_PRECISION == 2
    #define REAL(x) x
    #elif GEOGRAPHICLIB_PRECISION == 3
    #define REAL(x) x##L
    #elif GEOGRAPHICLIB_PRECISION == 4
    #define REAL(x) x##Q
    #else
    #define REAL(x) real(#x)
    #endif 

    and then use

    real f = 1/REAL(298.257223563); 

  • Perhaps use local static declarations to avoid the overhead of reevaluating constants, e.g.,

    static inline real pi() {
      using std::atan2;
      // pi is computed just once
      static const real pi = atan2(real(0), real(-1));
      return pi;
    } 

    This is not necessary for built-in floating point types, since the atan2 function will be evaluated at compile time.
  • In Utility::readarray and Utility::writearray, arrays of reals were treated as plain old data. This assumption now no longer holds and these functions needed special treatment.
  • volatile declarations don't apply.
• Phase 3, changes to support arbitrary precision which can be set at runtime.
  • The big change now is that the precision is not known at compile time. All static initializations which involve floating point numbers need to be eliminated.
    • Some static variables (e.g., tolerances which are expressed in terms of numeric_limits<double>::epsilon()) were made member variables and so initialized when the constructor was called.
    • Some simple static real arrays (e.g., the interpolating stencils for Geoid) were changed into integer arrays.
    • Some static variables where converted to static functions similar to the definition of pi() above.
    • All the static instances of classes where converted as follows

      // declaration
      static const Geodesic WGS84;
      // definition
      const Geodesic Geodesic::WGS84(Constants::WGS84_a(),
                                     Constants::WGS84_f());
      // use
      const Geodesic& geod = Geodesic::WGS84; 

      becomes

      // declaration
      static const Geodesic& WGS84();
      // definition
      const Geodesic& Geodesic::WGS84() {
        static const Geodesic wgs84(Constants::WGS84_a(),
                                    Constants::WGS84_f());
        return wgs84;
        static const Geodesic& WGS84();
      }
      // use
      const Geodesic& geod = Geodesic::WGS84(); 

      This is the so-called "construct on first use idiom". This is the most disruptive of the changes since it requires a different calling convention in user code. However the old static initializations were invoked every time a code linking to GeographicLib was started, even if the objects were not subsequently used. The new method only initializes the static objects if they are used.
  • numeric_limits<double>::digits is no longer a compile-time constant. It becomes numeric_limits<double>::digits().
  • Depending on the precision cos(π/2) might be negative. Similarly atan(tan(π/2)) may evaluate to −π/2. GeographicLib already handled this, because this happens with long doubles (64 bits in the fraction).
  • The precision needs to be set in each thread in a multi-processing applications (for an example, see examples/GeoidToGTX.cpp).
  • Passing numbers to functions by value incurs a substantially higher overhead than with doubles. This could be avoided by passing such arguments by reference. This was not done here because it would junk up the code to benefit a narrow application.
  • The constants in GeodesicExact, e.g., 831281402884796906843926125, can't be specified as long doubles nor long longs (since they have more than 64 significant bits):
    • first tried macro which expanded to a string (see the macro REAL above);
    • now use inline function to combine two long long ints each with at most 52 significant bits;
    • also needed to simplify one routine in GeodesicExact which took inordinately long (15 minutes) to compile using g++.

Back to Auxiliary latitudes. Forward to Change log. Up to Contents.