/usr/share/doc/geographiclib/html/highprec.html is in geographiclib-doc 1.49-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 | <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<meta name="generator" content="Doxygen 1.8.13"/>
<meta name="viewport" content="width=device-width, initial-scale=1"/>
<title>GeographicLib: Support for high precision arithmetic</title>
<link href="tabs.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="dynsections.js"></script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
extensions: ["tex2jax.js"],
jax: ["input/TeX","output/HTML-CSS"],
});
</script><script type="text/javascript" src="/usr/share/javascript/mathjax/MathJax.js/MathJax.js"></script>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
</head>
<body>
<div id="top"><!-- do not remove this div, it is closed by doxygen! -->
<div id="titlearea">
<table cellspacing="0" cellpadding="0">
<tbody>
<tr style="height: 56px;">
<td id="projectalign" style="padding-left: 0.5em;">
<div id="projectname">GeographicLib
 <span id="projectnumber">1.49</span>
</div>
</td>
</tr>
</tbody>
</table>
</div>
<!-- end header part -->
<!-- Generated by Doxygen 1.8.13 -->
<script type="text/javascript" src="menudata.js"></script>
<script type="text/javascript" src="menu.js"></script>
<script type="text/javascript">
$(function() {
initMenu('',false,false,'search.php','Search');
});
</script>
<div id="main-nav"></div>
</div><!-- top -->
<div class="header">
<div class="headertitle">
<div class="title">Support for high precision arithmetic </div> </div>
</div><!--header-->
<div class="contents">
<div class="textblock"><center> Back to <a class="el" href="auxlat.html">Auxiliary latitudes</a>. Forward to <a class="el" href="changes.html">Change log</a>. Up to <a class="el" href="index.html#contents">Contents</a>. </center><p>One of the goals with the algorithms in <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a> 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 <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a> by compiling it with <code>GEOGRAPHICLIB_PRECISION</code> equal to 4 or 5.</p>
<p>Here's what you should know:</p><ul>
<li>This is mainly for use for algorithm developers. It's not recommended for installation for all users on a system.</li>
<li>Configuring with <code>-D GEOGRAPHICLIB_PRECISION=4</code> gives quad precision (113-bit precision) via boost::multiprecision::float128; this requires:<ul>
<li><a href="http://www.boost.org">Boost</a>, version 1.54 or later,</li>
<li>the quadmath library (the package names are libquadmath and libquadmath-devel),</li>
<li>the use of g++.</li>
</ul>
</li>
<li>Configuring with <code>-D GEOGRAPHICLIB_PRECISION=5</code> gives arbitrary precision via mpfr::mpreal; this requires:<ul>
<li><a href="http://www.mpfr.org">MPFR</a>, version 3.0 or later,</li>
<li><a href="http://www.holoborodko.com/pavel/mpfr">MPFR C++</a> (version 3.6.5, dated 2016-12-19, or later),</li>
<li>a compiler which supports the explicit cast operator (e.g., g++ 4.5 or later, Visual Studio 12 2013 or later).</li>
</ul>
</li>
<li>MPFR, MPFR C++, and Boost all come with their own licenses. Be sure to respect these.</li>
<li>The indicated precision is used for <b>all</b> 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.</li>
<li>With MPFR, the precision should be set (using <a class="el" href="classGeographicLib_1_1Utility.html#ab95e26081d2c69019703adb8a4e4d3d0">Utility::set_digits</a>) just once before any other <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a> routines are called. Calling this function after other <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a> routines will lead to inconsistent results (because the precision of some constants like <a class="el" href="classGeographicLib_1_1Math.html#aca1580c771e7019bfe826512fba9b2f1">Math::pi()</a> is set when the functions are first called).</li>
<li>All the <a class="el" href="utilities.html">Utility programs</a> call <a class="el" href="classGeographicLib_1_1Utility.html#ab95e26081d2c69019703adb8a4e4d3d0">Utility::set_digits()</a> (with no arguments). This causes the precision (in bits) to be determined by the <code>GEOGRAPHICLIB_DIGITS</code> environment variable. If this is not defined the precision is set to 256 bits (about 77 decimal digits).</li>
<li>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:<ul>
<li>Round-off errors; these are reliably reduced when the precision is increased. For the most part, the algorithms used by <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a> are carefully tuned to minimize round-off errors, so that only a few bits of accuracy are lost.</li>
<li>Convergence errors due to not iterating certain algorithms to convergence. However, all iterative methods used by <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a> 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.</li>
<li>Truncation errors. Some classes (namely, <a class="el" href="classGeographicLib_1_1Geodesic.html" title="Geodesic calculations ">Geodesic</a> and <a class="el" href="classGeographicLib_1_1TransverseMercator.html" title="Transverse Mercator projection. ">TransverseMercator</a>) 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, <a class="el" href="classGeographicLib_1_1GeodesicExact.html" title="Exact geodesic calculations. ">GeodesicExact</a> and <a class="el" href="classGeographicLib_1_1TransverseMercatorExact.html" title="An exact implementation of the transverse Mercator projection. ">TransverseMercatorExact</a> are somewhat slower classes offering the same functionality implemented with <a class="el" href="classGeographicLib_1_1EllipticFunction.html" title="Elliptic integrals and functions. ">EllipticFunction</a>. These classes provide arbitrary accuracy. (However, a caveat is that the evaluation of the area in <a class="el" href="classGeographicLib_1_1GeodesicExact.html" title="Exact geodesic calculations. ">GeodesicExact</a> still uses a series (albeit of considerably higher order). So the area calculations are always have a finite truncation error.)</li>
<li>Quantization errors. <a class="el" href="classGeographicLib_1_1Geoid.html" title="Looking up the height of the geoid above the ellipsoid. ">Geoid</a>, <a class="el" href="classGeographicLib_1_1GravityModel.html" title="Model of the earth's gravity field. ">GravityModel</a>, and <a class="el" href="classGeographicLib_1_1MagneticModel.html" title="Model of the earth's magnetic field. ">MagneticModel</a> all depend on external data files. The coefficient files for <a class="el" href="classGeographicLib_1_1GravityModel.html" title="Model of the earth's gravity field. ">GravityModel</a> and <a class="el" href="classGeographicLib_1_1MagneticModel.html" title="Model of the earth's magnetic field. ">MagneticModel</a> store the coefficients as IEEE doubles (and perhaps these coefficients can be regarded as exact). However, with <a class="el" href="classGeographicLib_1_1Geoid.html" title="Looking up the height of the geoid above the ellipsoid. ">Geoid</a>, 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 <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a>, e.g., the flattening of the WGS84 ellipsoid, are evaluated as exact decimal numbers.</li>
</ul>
</li>
<li>Where might high accuracy be important?<ul>
<li>checking the truncation error of series approximations;</li>
<li>checking for excessive round-off errors (typically due to subtraction);</li>
<li>checking the round-off error in computing areas of many-sided polygons;</li>
<li>checking the summation of high order spherical harmonic expansions (where underflow and overflow may also be a problem).</li>
</ul>
</li>
<li>Because only a tiny number of people will be interested in using this facility:<ul>
<li>the cmake support for the required libraries is rudimentary;</li>
<li>however geographiclib-config.cmake does export <code>GEOGRAPHICLIB_PRECISION</code> and <code>GEOGRAPHICLIB_HIGHPREC_LIBRARIES</code>, the libraries providing the support for high-precision arithmetic;</li>
<li>support for the C++11 mathematical functions and the explicit cast operator is required;</li>
<li>quad precision is only available on Linux;</li>
<li>mpfr has been mostly tested on Linux (but it works on Windows with Visual Studio 12 and MacOS too).</li>
</ul>
</li>
</ul>
<p>The following steps needed to be taken</p>
<ul>
<li>Phase 1, make sure you can switch easily between double, float, and long double.<ul>
<li>use <code>#include <cmath></code> instead of <code>#include <math.h></code>;</li>
<li>use, e.g., <code>std::sqrt</code> instead of <code>sqrt</code> in header files (similarly for <code>sin</code>, <code>cos</code>, <code>atan2</code>, etc.);</li>
<li>use <code>using namespace std;</code> and plain <code>sqrt</code>, etc., in code files;</li>
<li>express all convergence criteria in terms of<div class="fragment"><div class="line">numeric_limits<double>::epsilon() </div></div><!-- fragment --> etc., instead of using "magic constants", such as 1.0e-15;</li>
<li>use <code>typedef double real;</code> and replace all occurrences of <code>double</code> by <code>real</code>;</li>
<li>write all literals by, e.g., <code>real(0.5)</code>. Some constants might need the L suffix, e.g., <code>real f = 1/real(298.257223563L)</code> (but see below);</li>
<li>Change the typedef of <code>real</code> to <code>float</code> or <code>long double</code>, compile, and test. In this way, the library can be run with any of the three basic floating point types.</li>
<li>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 <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a> 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.</li>
</ul>
</li>
<li>Phase 2, changes to support arbitrary, but fixed, precision<ul>
<li>Use, e.g.,<div class="fragment"><div class="line"><span class="keyword">typedef</span> boost::multiprecision::float128 <a class="code" href="GeodSolve_8cpp.html#a5caf95d46b184d9ca1d3764b3781b3c9">real</a>; </div></div><!-- fragment --></li>
<li>Change <code>std::sqrt(...)</code>, etc. to<div class="fragment"><div class="line"><span class="keyword">using</span> std::sqrt;</div><div class="line">sqrt(...) </div></div><!-- fragment --> (but note that <code>std::max</code> can stay). It's only necessary to do this in header files (code files already have <code>using namespace std;</code>).</li>
<li>In the case of boost's multiprecision numbers, the C++11 mathematical functions need special treatment, see <a class="el" href="Math_8hpp.html" title="Header for GeographicLib::Math class. ">Math.hpp</a>.</li>
<li>If necessary, use series with additional terms to improve the accuracy.</li>
<li>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 <a class="el" href="classGeographicLib_1_1NormalGravity.html" title="The normal gravity of the earth. ">NormalGravity</a> converged too slowly; this was replaced by Newton's method.</li>
<li>If necessary, increase the maximum allowed iteration count in root finding loops. Also throw an exception of the maximum iteration count is exceeded.</li>
<li>Write literal constants in a way that works for any precision, e.g., <div class="fragment"><div class="line">real f = 1/( <a class="code" href="GeodSolve_8cpp.html#a5caf95d46b184d9ca1d3764b3781b3c9">real</a>(298257223563LL) / 1000000000 ); </div></div><!-- fragment --> [Note that<div class="fragment"><div class="line">real f = 1/( 298 + <a class="code" href="GeodSolve_8cpp.html#a5caf95d46b184d9ca1d3764b3781b3c9">real</a>(257223563) / 1000000000 ); </div></div><!-- fragment --> and <code>1/real(298.257223563L)</code> are susceptible to double rounding errors. We normally want to avoid such errors when real is a double.]</li>
<li>For arbitrary constants, you might have to resort to macros<div class="fragment"><div class="line"><span class="preprocessor">#if GEOGRAPHICLIB_PRECISION == 1</span></div><div class="line"><span class="preprocessor">#define REAL(x) x##F</span></div><div class="line"><span class="preprocessor">#elif GEOGRAPHICLIB_PRECISION == 2</span></div><div class="line"><span class="preprocessor">#define REAL(x) x</span></div><div class="line"><span class="preprocessor">#elif GEOGRAPHICLIB_PRECISION == 3</span></div><div class="line"><span class="preprocessor">#define REAL(x) x##L</span></div><div class="line"><span class="preprocessor">#elif GEOGRAPHICLIB_PRECISION == 4</span></div><div class="line"><span class="preprocessor">#define REAL(x) x##Q</span></div><div class="line"><span class="preprocessor">#else</span></div><div class="line"><span class="preprocessor">#define REAL(x) real(#x)</span></div><div class="line"><span class="preprocessor">#endif </span></div></div><!-- fragment --> and then use<div class="fragment"><div class="line">real f = 1/REAL(298.257223563); </div></div><!-- fragment --></li>
<li>Perhaps use local static declarations to avoid the overhead of reevaluating constants, e.g.,<div class="fragment"><div class="line"><span class="keyword">static</span> <span class="keyword">inline</span> real pi() {</div><div class="line"> <span class="keyword">using</span> std::atan2;</div><div class="line"> <span class="comment">// pi is computed just once</span></div><div class="line"> <span class="keyword">static</span> <span class="keyword">const</span> real pi = atan2(<a class="code" href="GeodSolve_8cpp.html#a5caf95d46b184d9ca1d3764b3781b3c9">real</a>(0), <a class="code" href="GeodSolve_8cpp.html#a5caf95d46b184d9ca1d3764b3781b3c9">real</a>(-1));</div><div class="line"> <span class="keywordflow">return</span> pi;</div><div class="line">} </div></div><!-- fragment --> This is <b>not</b> necessary for built-in floating point types, since the atan2 function will be evaluated at compile time.</li>
<li>In <a class="el" href="classGeographicLib_1_1Utility.html#a23e88040ceb60bd3fe28efc04f9119f8">Utility::readarray</a> and <a class="el" href="classGeographicLib_1_1Utility.html#a276eb20ace49c5260e1878c7d8aefd31">Utility::writearray</a>, arrays of reals were treated as plain old data. This assumption now no longer holds and these functions needed special treatment.</li>
<li>volatile declarations don't apply.</li>
</ul>
</li>
<li>Phase 3, changes to support arbitrary precision which can be set at runtime.<ul>
<li>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.<ul>
<li>Some static variables (e.g., tolerances which are expressed in terms of <code>numeric_limits<double>::epsilon()</code>) were made member variables and so initialized when the constructor was called.</li>
<li>Some simple static real arrays (e.g., the interpolating stencils for <a class="el" href="classGeographicLib_1_1Geoid.html" title="Looking up the height of the geoid above the ellipsoid. ">Geoid</a>) were changed into integer arrays.</li>
<li>Some static variables where converted to static functions similar to the definition of pi() above.</li>
<li>All the static instances of classes where converted as follows <div class="fragment"><div class="line"><span class="comment">// declaration</span></div><div class="line"><span class="keyword">static</span> <span class="keyword">const</span> Geodesic WGS84;</div><div class="line"><span class="comment">// definition</span></div><div class="line"><span class="keyword">const</span> Geodesic <a class="code" href="classGeographicLib_1_1Geodesic.html#a479d30b83e1d2d32ad31c96b4df1543a">Geodesic::WGS84</a>(<a class="code" href="classGeographicLib_1_1Constants.html#ae12127984ac6713823746d917b4abfa7">Constants::WGS84_a</a>(),</div><div class="line"> <a class="code" href="classGeographicLib_1_1Constants.html#acc5244425bb484594be51b27d56cd335">Constants::WGS84_f</a>());</div><div class="line"><span class="comment">// use</span></div><div class="line"><span class="keyword">const</span> Geodesic& geod = <a class="code" href="classGeographicLib_1_1Geodesic.html#a479d30b83e1d2d32ad31c96b4df1543a">Geodesic::WGS84</a>; </div></div><!-- fragment --> becomes<div class="fragment"><div class="line"><span class="comment">// declaration</span></div><div class="line"><span class="keyword">static</span> <span class="keyword">const</span> Geodesic& WGS84();</div><div class="line"><span class="comment">// definition</span></div><div class="line"><span class="keyword">const</span> Geodesic& <a class="code" href="classGeographicLib_1_1Geodesic.html#a479d30b83e1d2d32ad31c96b4df1543a">Geodesic::WGS84</a>() {</div><div class="line"> <span class="keyword">static</span> <span class="keyword">const</span> <a class="code" href="classGeographicLib_1_1Geodesic.html#a455300c36e6caa70968115416e1573a4">Geodesic</a> wgs84(<a class="code" href="classGeographicLib_1_1Constants.html#ae12127984ac6713823746d917b4abfa7">Constants::WGS84_a</a>(),</div><div class="line"> <a class="code" href="classGeographicLib_1_1Constants.html#acc5244425bb484594be51b27d56cd335">Constants::WGS84_f</a>());</div><div class="line"> <span class="keywordflow">return</span> wgs84;</div><div class="line"> <span class="keyword">static</span> <span class="keyword">const</span> <a class="code" href="classGeographicLib_1_1Geodesic.html#a455300c36e6caa70968115416e1573a4">Geodesic</a>& <a class="code" href="classGeographicLib_1_1Geodesic.html#a479d30b83e1d2d32ad31c96b4df1543a">WGS84</a>();</div><div class="line">}</div><div class="line"><span class="comment">// use</span></div><div class="line"><span class="keyword">const</span> <a class="code" href="classGeographicLib_1_1Geodesic.html#a455300c36e6caa70968115416e1573a4">Geodesic</a>& geod = <a class="code" href="classGeographicLib_1_1Geodesic.html#a479d30b83e1d2d32ad31c96b4df1543a">Geodesic::WGS84</a>(); </div></div><!-- fragment --> This is the so-called <a href="https://isocpp.org/wiki/faq/ctors#static-init-order-on-first-use">"construct on first use idiom"</a>. 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 <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a> was started, even if the objects were not subsequently used. The new method only initializes the static objects if they are used.</li>
</ul>
</li>
<li><code>numeric_limits<double>::digits</code> is no longer a compile-time constant. It becomes <code>numeric_limits<double>::digits()</code>.</li>
<li>Depending on the precision cos(π/2) might be negative. Similarly atan(tan(π/2)) may evaluate to −π/2. <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a> already handled this, because this happens with long doubles (64 bits in the fraction).</li>
<li>The precision needs to be set in each thread in a multi-processing applications (for an example, see <code>examples/GeoidToGTX.cpp</code>).</li>
<li>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 <b>not</b> done here because it would junk up the code to benefit a narrow application.</li>
<li>The constants in <a class="el" href="classGeographicLib_1_1GeodesicExact.html" title="Exact geodesic calculations. ">GeodesicExact</a>, e.g., 831281402884796906843926125, can't be specified as long doubles nor long longs (since they have more than 64 significant bits):<ul>
<li>first tried macro which expanded to a string (see the macro REAL above);</li>
<li>now use inline function to combine two long long ints each with at most 52 significant bits;</li>
<li>also needed to simplify one routine in <a class="el" href="classGeographicLib_1_1GeodesicExact.html" title="Exact geodesic calculations. ">GeodesicExact</a> which took inordinately long (15 minutes) to compile using g++;</li>
<li>but Visual Studio 12 then complained (with an <a href="http://connect.microsoft.com/VisualStudio/feedback/details/920594">internal compiler error</a> presumably due to overly aggressive whole-file optimization); fixed by splitting one function into a separate file.</li>
</ul>
</li>
</ul>
</li>
</ul>
<center> Back to <a class="el" href="auxlat.html">Auxiliary latitudes</a>. Forward to <a class="el" href="changes.html">Change log</a>. Up to <a class="el" href="index.html#contents">Contents</a>. </center> </div></div><!-- contents -->
<!-- start footer part -->
<hr class="footer"/><address class="footer"><small>
Generated by  <a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/>
</a> 1.8.13
</small></address>
</body>
</html>
|