GeographicLib  2.1
Math.hpp
Go to the documentation of this file.
1 /**
2  * \file Math.hpp
3  * \brief Header for GeographicLib::Math class
4  *
5  * Copyright (c) Charles Karney (2008-2022) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
10 // Constants.hpp includes Math.hpp. Place this include outside Math.hpp's
11 // include guard to enforce this ordering.
13 
14 #if !defined(GEOGRAPHICLIB_MATH_HPP)
15 #define GEOGRAPHICLIB_MATH_HPP 1
16 
17 #if !defined(GEOGRAPHICLIB_WORDS_BIGENDIAN)
18 # define GEOGRAPHICLIB_WORDS_BIGENDIAN 0
19 #endif
20 
21 #if !defined(GEOGRAPHICLIB_HAVE_LONG_DOUBLE)
22 # define GEOGRAPHICLIB_HAVE_LONG_DOUBLE 0
23 #endif
24 
25 #if !defined(GEOGRAPHICLIB_PRECISION)
26 /**
27  * The precision of floating point numbers used in %GeographicLib. 1 means
28  * float (single precision); 2 (the default) means double; 3 means long double;
29  * 4 is reserved for quadruple precision. Nearly all the testing has been
30  * carried out with doubles and that's the recommended configuration. In order
31  * for long double to be used, GEOGRAPHICLIB_HAVE_LONG_DOUBLE needs to be
32  * defined. Note that with Microsoft Visual Studio, long double is the same as
33  * double.
34  **********************************************************************/
35 # define GEOGRAPHICLIB_PRECISION 2
36 #endif
37 
38 #include <cmath>
39 #include <algorithm>
40 #include <limits>
41 
42 #if GEOGRAPHICLIB_PRECISION == 4
43 #include <boost/version.hpp>
44 #include <boost/multiprecision/float128.hpp>
45 #include <boost/math/special_functions.hpp>
46 #elif GEOGRAPHICLIB_PRECISION == 5
47 #include <mpreal.h>
48 #endif
49 
50 #if GEOGRAPHICLIB_PRECISION > 3
51 // volatile keyword makes no sense for multiprec types
52 #define GEOGRAPHICLIB_VOLATILE
53 // Signal a convergence failure with multiprec types by throwing an exception
54 // at loop exit.
55 #define GEOGRAPHICLIB_PANIC \
56  (throw GeographicLib::GeographicErr("Convergence failure"), false)
57 #else
58 #define GEOGRAPHICLIB_VOLATILE volatile
59 // Ignore convergence failures with standard floating points types by allowing
60 // loop to exit cleanly.
61 #define GEOGRAPHICLIB_PANIC false
62 #endif
63 
64 namespace GeographicLib {
65 
66  /**
67  * \brief Mathematical functions needed by %GeographicLib
68  *
69  * Define mathematical functions in order to localize system dependencies and
70  * to provide generic versions of the functions. In addition define a real
71  * type to be used by %GeographicLib.
72  *
73  * Example of use:
74  * \include example-Math.cpp
75  **********************************************************************/
77  private:
78  void dummy(); // Static check for GEOGRAPHICLIB_PRECISION
79  Math() = delete; // Disable constructor
80  public:
81 
82 #if GEOGRAPHICLIB_HAVE_LONG_DOUBLE
83  /**
84  * The extended precision type for real numbers, used for some testing.
85  * This is long double on computers with this type; otherwise it is double.
86  **********************************************************************/
87  typedef long double extended;
88 #else
89  typedef double extended;
90 #endif
91 
92 #if GEOGRAPHICLIB_PRECISION == 2
93  /**
94  * The real type for %GeographicLib. Nearly all the testing has been done
95  * with \e real = double. However, the algorithms should also work with
96  * float and long double (where available). (<b>CAUTION</b>: reasonable
97  * accuracy typically cannot be obtained using floats.)
98  **********************************************************************/
99  typedef double real;
100 #elif GEOGRAPHICLIB_PRECISION == 1
101  typedef float real;
102 #elif GEOGRAPHICLIB_PRECISION == 3
103  typedef extended real;
104 #elif GEOGRAPHICLIB_PRECISION == 4
105  typedef boost::multiprecision::float128 real;
106 #elif GEOGRAPHICLIB_PRECISION == 5
107  typedef mpfr::mpreal real;
108 #else
109  typedef double real;
110 #endif
111 
112  /**
113  * The constants defining the standard (Babylonian) meanings of degrees,
114  * minutes, and seconds, for angles. Read the constants as follows (for
115  * example): \e ms = 60 is the ratio 1 minute / 1 second. The
116  * abbreviations are
117  * - \e t a whole turn (360&deg;)
118  * - \e h a half turn (180&deg;)
119  * - \e q a quarter turn (a right angle = 90&deg;)
120  * - \e d a degree
121  * - \e m a minute
122  * - \e s a second
123  * .
124  * Note that degree() is ratio 1 degree / 1 radian, thus, for example,
125  * Math::degree() * Math::qd is the ratio 1 quarter turn / 1 radian =
126  * &pi;/2.
127  *
128  * Defining all these in one place would mean that it's simple to convert
129  * to the centesimal system for measuring angles. The DMS class assumes
130  * that Math::dm and Math::ms are less than or equal to 100 (so that two
131  * digits suffice for the integer parts of the minutes and degrees
132  * components of an angle). Switching to the centesimal convention will
133  * break most of the tests. Also the normal definition of degree is baked
134  * into some classes, e.g., UTMUPS, MGRS, Georef, Geohash, etc.
135  **********************************************************************/
136 #if GEOGRAPHICLIB_PRECISION == 4
137  static const int
138 #else
139  enum dms {
140 #endif
141  qd = 90, ///< degrees per quarter turn
142  dm = 60, ///< minutes per degree
143  ms = 60, ///< seconds per minute
144  hd = 2 * qd, ///< degrees per half turn
145  td = 2 * hd, ///< degrees per turn
146  ds = dm * ms ///< seconds per degree
147 #if GEOGRAPHICLIB_PRECISION == 4
148  ;
149 #else
150  };
151 #endif
152 
153  /**
154  * @return the number of bits of precision in a real number.
155  **********************************************************************/
156  static int digits();
157 
158  /**
159  * Set the binary precision of a real number.
160  *
161  * @param[in] ndigits the number of bits of precision.
162  * @return the resulting number of bits of precision.
163  *
164  * This only has an effect when GEOGRAPHICLIB_PRECISION = 5. See also
165  * Utility::set_digits for caveats about when this routine should be
166  * called.
167  **********************************************************************/
168  static int set_digits(int ndigits);
169 
170  /**
171  * @return the number of decimal digits of precision in a real number.
172  **********************************************************************/
173  static int digits10();
174 
175  /**
176  * Number of additional decimal digits of precision for real relative to
177  * double (0 for float).
178  **********************************************************************/
179  static int extra_digits();
180 
181  /**
182  * true if the machine is big-endian.
183  **********************************************************************/
184  static const bool bigendian = GEOGRAPHICLIB_WORDS_BIGENDIAN;
185 
186  /**
187  * @tparam T the type of the returned value.
188  * @return &pi;.
189  **********************************************************************/
190  template<typename T = real> static T pi() {
191  using std::atan2;
192  static const T pi = atan2(T(0), T(-1));
193  return pi;
194  }
195 
196  /**
197  * @tparam T the type of the returned value.
198  * @return the number of radians in a degree.
199  **********************************************************************/
200  template<typename T = real> static T degree() {
201  static const T degree = pi<T>() / T(hd);
202  return degree;
203  }
204 
205  /**
206  * Square a number.
207  *
208  * @tparam T the type of the argument and the returned value.
209  * @param[in] x
210  * @return <i>x</i><sup>2</sup>.
211  **********************************************************************/
212  template<typename T> static T sq(T x)
213  { return x * x; }
214 
215  /**
216  * Normalize a two-vector.
217  *
218  * @tparam T the type of the argument and the returned value.
219  * @param[in,out] x on output set to <i>x</i>/hypot(<i>x</i>, <i>y</i>).
220  * @param[in,out] y on output set to <i>y</i>/hypot(<i>x</i>, <i>y</i>).
221  **********************************************************************/
222  template<typename T> static void norm(T& x, T& y) {
223 #if defined(_MSC_VER) && defined(_M_IX86)
224  // hypot for Visual Studio (A=win32) fails monotonicity, e.g., with
225  // x = 0.6102683302836215
226  // y1 = 0.7906090004346522
227  // y2 = y1 + 1e-16
228  // the test
229  // hypot(x, y2) >= hypot(x, y1)
230  // fails. Reported 2021-03-14:
231  // https://developercommunity.visualstudio.com/t/1369259
232  // See also:
233  // https://bugs.python.org/issue43088
234  using std::sqrt; T h = sqrt(x * x + y * y);
235 #else
236  using std::hypot; T h = hypot(x, y);
237 #endif
238  x /= h; y /= h;
239  }
240 
241  /**
242  * The error-free sum of two numbers.
243  *
244  * @tparam T the type of the argument and the returned value.
245  * @param[in] u
246  * @param[in] v
247  * @param[out] t the exact error given by (\e u + \e v) - \e s.
248  * @return \e s = round(\e u + \e v).
249  *
250  * See D. E. Knuth, TAOCP, Vol 2, 4.2.2, Theorem B.
251  *
252  * \note \e t can be the same as one of the first two arguments.
253  **********************************************************************/
254  template<typename T> static T sum(T u, T v, T& t);
255 
256  /**
257  * Evaluate a polynomial.
258  *
259  * @tparam T the type of the arguments and returned value.
260  * @param[in] N the order of the polynomial.
261  * @param[in] p the coefficient array (of size \e N + 1) with
262  * <i>p</i><sub>0</sub> being coefficient of <i>x</i><sup><i>N</i></sup>.
263  * @param[in] x the variable.
264  * @return the value of the polynomial.
265  *
266  * Evaluate &sum;<sub><i>n</i>=0..<i>N</i></sub>
267  * <i>p</i><sub><i>n</i></sub> <i>x</i><sup><i>N</i>&minus;<i>n</i></sup>.
268  * Return 0 if \e N &lt; 0. Return <i>p</i><sub>0</sub>, if \e N = 0 (even
269  * if \e x is infinite or a nan). The evaluation uses Horner's method.
270  **********************************************************************/
271  template<typename T> static T polyval(int N, const T p[], T x) {
272  // This used to employ Math::fma; but that's too slow and it seemed not
273  // to improve the accuracy noticeably. This might change when there's
274  // direct hardware support for fma.
275  T y = N < 0 ? 0 : *p++;
276  while (--N >= 0) y = y * x + *p++;
277  return y;
278  }
279 
280  /**
281  * Normalize an angle.
282  *
283  * @tparam T the type of the argument and returned value.
284  * @param[in] x the angle in degrees.
285  * @return the angle reduced to the range [&minus;180&deg;, 180&deg;].
286  *
287  * The range of \e x is unrestricted. If the result is &plusmn;0&deg; or
288  * &plusmn;180&deg; then the sign is the sign of \e x.
289  **********************************************************************/
290  template<typename T> static T AngNormalize(T x);
291 
292  /**
293  * Normalize a latitude.
294  *
295  * @tparam T the type of the argument and returned value.
296  * @param[in] x the angle in degrees.
297  * @return x if it is in the range [&minus;90&deg;, 90&deg;], otherwise
298  * return NaN.
299  **********************************************************************/
300  template<typename T> static T LatFix(T x)
301  { using std::fabs; return fabs(x) > T(qd) ? NaN<T>() : x; }
302 
303  /**
304  * The exact difference of two angles reduced to
305  * [&minus;180&deg;, 180&deg;].
306  *
307  * @tparam T the type of the arguments and returned value.
308  * @param[in] x the first angle in degrees.
309  * @param[in] y the second angle in degrees.
310  * @param[out] e the error term in degrees.
311  * @return \e d, the truncated value of \e y &minus; \e x.
312  *
313  * This computes \e z = \e y &minus; \e x exactly, reduced to
314  * [&minus;180&deg;, 180&deg;]; and then sets \e z = \e d + \e e where \e d
315  * is the nearest representable number to \e z and \e e is the truncation
316  * error. If \e z = &plusmn;0&deg; or &plusmn;180&deg;, then the sign of
317  * \e d is given by the sign of \e y &minus; \e x. The maximum absolute
318  * value of \e e is 2<sup>&minus;26</sup> (for doubles).
319  **********************************************************************/
320  template<typename T> static T AngDiff(T x, T y, T& e);
321 
322  /**
323  * Difference of two angles reduced to [&minus;180&deg;, 180&deg;]
324  *
325  * @tparam T the type of the arguments and returned value.
326  * @param[in] x the first angle in degrees.
327  * @param[in] y the second angle in degrees.
328  * @return \e y &minus; \e x, reduced to the range [&minus;180&deg;,
329  * 180&deg;].
330  *
331  * The result is equivalent to computing the difference exactly, reducing
332  * it to [&minus;180&deg;, 180&deg;] and rounding the result.
333  **********************************************************************/
334  template<typename T> static T AngDiff(T x, T y)
335  { T e; return AngDiff(x, y, e); }
336 
337  /**
338  * Coarsen a value close to zero.
339  *
340  * @tparam T the type of the argument and returned value.
341  * @param[in] x
342  * @return the coarsened value.
343  *
344  * The makes the smallest gap in \e x = 1/16 &minus; nextafter(1/16, 0) =
345  * 1/2<sup>57</sup> for doubles = 0.8 pm on the earth if \e x is an angle
346  * in degrees. (This is about 2000 times more resolution than we get with
347  * angles around 90&deg;.) We use this to avoid having to deal with near
348  * singular cases when \e x is non-zero but tiny (e.g.,
349  * 10<sup>&minus;200</sup>). This sign of &plusmn;0 is preserved.
350  **********************************************************************/
351  template<typename T> static T AngRound(T x);
352 
353  /**
354  * Evaluate the sine and cosine function with the argument in degrees
355  *
356  * @tparam T the type of the arguments.
357  * @param[in] x in degrees.
358  * @param[out] sinx sin(<i>x</i>).
359  * @param[out] cosx cos(<i>x</i>).
360  *
361  * The results obey exactly the elementary properties of the trigonometric
362  * functions, e.g., sin 9&deg; = cos 81&deg; = &minus; sin 123456789&deg;.
363  * If x = &minus;0 or a negative multiple of 180&deg;, then \e sinx =
364  * &minus;0; this is the only case where &minus;0 is returned.
365  **********************************************************************/
366  template<typename T> static void sincosd(T x, T& sinx, T& cosx);
367 
368  /**
369  * Evaluate the sine and cosine with reduced argument plus correction
370  *
371  * @tparam T the type of the arguments.
372  * @param[in] x reduced angle in degrees.
373  * @param[in] t correction in degrees.
374  * @param[out] sinx sin(<i>x</i> + <i>t</i>).
375  * @param[out] cosx cos(<i>x</i> + <i>t</i>).
376  *
377  * This is a variant of Math::sincosd allowing a correction to the angle to
378  * be supplied. \e x must be in [&minus;180&deg;, 180&deg;] and \e t is
379  * assumed to be a <i>small</i> correction. Math::AngRound is applied to
380  * the reduced angle to prevent problems with \e x + \e t being extremely
381  * close but not exactly equal to one of the four cardinal directions.
382  **********************************************************************/
383  template<typename T> static void sincosde(T x, T t, T& sinx, T& cosx);
384 
385  /**
386  * Evaluate the sine function with the argument in degrees
387  *
388  * @tparam T the type of the argument and the returned value.
389  * @param[in] x in degrees.
390  * @return sin(<i>x</i>).
391  *
392  * The result is +0 for \e x = +0 and positive multiples of 180&deg;. The
393  * result is &minus;0 for \e x = -0 and negative multiples of 180&deg;.
394  **********************************************************************/
395  template<typename T> static T sind(T x);
396 
397  /**
398  * Evaluate the cosine function with the argument in degrees
399  *
400  * @tparam T the type of the argument and the returned value.
401  * @param[in] x in degrees.
402  * @return cos(<i>x</i>).
403  *
404  * The result is +0 for \e x an odd multiple of 90&deg;.
405  **********************************************************************/
406  template<typename T> static T cosd(T x);
407 
408  /**
409  * Evaluate the tangent function with the argument in degrees
410  *
411  * @tparam T the type of the argument and the returned value.
412  * @param[in] x in degrees.
413  * @return tan(<i>x</i>).
414  *
415  * If \e x is an odd multiple of 90&deg;, then a suitably large (but
416  * finite) value is returned.
417  **********************************************************************/
418  template<typename T> static T tand(T x);
419 
420  /**
421  * Evaluate the atan2 function with the result in degrees
422  *
423  * @tparam T the type of the arguments and the returned value.
424  * @param[in] y
425  * @param[in] x
426  * @return atan2(<i>y</i>, <i>x</i>) in degrees.
427  *
428  * The result is in the range [&minus;180&deg; 180&deg;]. N.B.,
429  * atan2d(&plusmn;0, &minus;1) = &plusmn;180&deg;.
430  **********************************************************************/
431  template<typename T> static T atan2d(T y, T x);
432 
433  /**
434  * Evaluate the atan function with the result in degrees
435  *
436  * @tparam T the type of the argument and the returned value.
437  * @param[in] x
438  * @return atan(<i>x</i>) in degrees.
439  **********************************************************************/
440  template<typename T> static T atand(T x);
441 
442  /**
443  * Evaluate <i>e</i> atanh(<i>e x</i>)
444  *
445  * @tparam T the type of the argument and the returned value.
446  * @param[in] x
447  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
448  * sqrt(|<i>e</i><sup>2</sup>|)
449  * @return <i>e</i> atanh(<i>e x</i>)
450  *
451  * If <i>e</i><sup>2</sup> is negative (<i>e</i> is imaginary), the
452  * expression is evaluated in terms of atan.
453  **********************************************************************/
454  template<typename T> static T eatanhe(T x, T es);
455 
456  /**
457  * tan&chi; in terms of tan&phi;
458  *
459  * @tparam T the type of the argument and the returned value.
460  * @param[in] tau &tau; = tan&phi;
461  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
462  * sqrt(|<i>e</i><sup>2</sup>|)
463  * @return &tau;&prime; = tan&chi;
464  *
465  * See Eqs. (7--9) of
466  * C. F. F. Karney,
467  * <a href="https://doi.org/10.1007/s00190-011-0445-3">
468  * Transverse Mercator with an accuracy of a few nanometers,</a>
469  * J. Geodesy 85(8), 475--485 (Aug. 2011)
470  * (preprint
471  * <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>).
472  **********************************************************************/
473  template<typename T> static T taupf(T tau, T es);
474 
475  /**
476  * tan&phi; in terms of tan&chi;
477  *
478  * @tparam T the type of the argument and the returned value.
479  * @param[in] taup &tau;&prime; = tan&chi;
480  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
481  * sqrt(|<i>e</i><sup>2</sup>|)
482  * @return &tau; = tan&phi;
483  *
484  * See Eqs. (19--21) of
485  * C. F. F. Karney,
486  * <a href="https://doi.org/10.1007/s00190-011-0445-3">
487  * Transverse Mercator with an accuracy of a few nanometers,</a>
488  * J. Geodesy 85(8), 475--485 (Aug. 2011)
489  * (preprint
490  * <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>).
491  **********************************************************************/
492  template<typename T> static T tauf(T taup, T es);
493 
494  /**
495  * The NaN (not a number)
496  *
497  * @tparam T the type of the returned value.
498  * @return NaN if available, otherwise return the max real of type T.
499  **********************************************************************/
500  template<typename T = real> static T NaN();
501 
502  /**
503  * Infinity
504  *
505  * @tparam T the type of the returned value.
506  * @return infinity if available, otherwise return the max real.
507  **********************************************************************/
508  template<typename T = real> static T infinity();
509 
510  /**
511  * Swap the bytes of a quantity
512  *
513  * @tparam T the type of the argument and the returned value.
514  * @param[in] x
515  * @return x with its bytes swapped.
516  **********************************************************************/
517  template<typename T> static T swab(T x) {
518  union {
519  T r;
520  unsigned char c[sizeof(T)];
521  } b;
522  b.r = x;
523  for (int i = sizeof(T)/2; i--; )
524  std::swap(b.c[i], b.c[sizeof(T) - 1 - i]);
525  return b.r;
526  }
527 
528  };
529 
530 } // namespace GeographicLib
531 
532 #endif // GEOGRAPHICLIB_MATH_HPP
Header for GeographicLib::Constants class.
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:67
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
#define GEOGRAPHICLIB_WORDS_BIGENDIAN
Definition: Math.hpp:18
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:76
static T degree()
Definition: Math.hpp:200
static T LatFix(T x)
Definition: Math.hpp:300
double extended
Definition: Math.hpp:89
static void norm(T &x, T &y)
Definition: Math.hpp:222
static T sq(T x)
Definition: Math.hpp:212
static T pi()
Definition: Math.hpp:190
static T polyval(int N, const T p[], T x)
Definition: Math.hpp:271
static T AngDiff(T x, T y)
Definition: Math.hpp:334
static T swab(T x)
Definition: Math.hpp:517
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)