GeographicLib  1.43
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-2015) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
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 /**
18  * Are C++11 math functions available?
19  **********************************************************************/
20 #if !defined(GEOGRAPHICLIB_CXX11_MATH)
21 // Recent versions of g++ -std=c++11 (4.7 and later?) set __cplusplus to 201103
22 // and support the new C++11 mathematical functions, std::atanh, etc. However
23 // the Android toolchain, which uses g++ -std=c++11 (4.8 as of 2014-03-11,
24 // according to Pullan Lu), does not support std::atanh. Android toolchains
25 // might define __ANDROID__ or ANDROID; so need to check both. With OSX the
26 // version is GNUC version 4.2 and __cplusplus is set to 201103, so remove the
27 // version check on GNUC.
28 # if defined(__GNUC__) && __cplusplus >= 201103 && \
29  !(defined(__ANDROID__) || defined(ANDROID) || defined(__CYGWIN__))
30 # define GEOGRAPHICLIB_CXX11_MATH 1
31 // Visual C++ 12 supports these functions
32 # elif defined(_MSC_VER) && _MSC_VER >= 1800
33 # define GEOGRAPHICLIB_CXX11_MATH 1
34 # else
35 # define GEOGRAPHICLIB_CXX11_MATH 0
36 # endif
37 #endif
38 
39 #if !defined(GEOGRAPHICLIB_WORDS_BIGENDIAN)
40 # define GEOGRAPHICLIB_WORDS_BIGENDIAN 0
41 #endif
42 
43 #if !defined(GEOGRAPHICLIB_HAVE_LONG_DOUBLE)
44 # define GEOGRAPHICLIB_HAVE_LONG_DOUBLE 0
45 #endif
46 
47 #if !defined(GEOGRAPHICLIB_PRECISION)
48 /**
49  * The precision of floating point numbers used in %GeographicLib. 1 means
50  * float (single precision); 2 (the default) means double; 3 means long double;
51  * 4 is reserved for quadruple precision. Nearly all the testing has been
52  * carried out with doubles and that's the recommended configuration. In order
53  * for long double to be used, GEOGRAPHICLIB_HAVE_LONG_DOUBLE needs to be
54  * defined. Note that with Microsoft Visual Studio, long double is the same as
55  * double.
56  **********************************************************************/
57 # define GEOGRAPHICLIB_PRECISION 2
58 #endif
59 
60 #include <cmath>
61 #include <algorithm>
62 #include <limits>
63 
64 #if GEOGRAPHICLIB_PRECISION == 4
65 #include <boost/version.hpp>
66 #if BOOST_VERSION >= 105600
67 #include <boost/cstdfloat.hpp>
68 #endif
69 #include <boost/multiprecision/float128.hpp>
70 #include <boost/math/special_functions.hpp>
71 __float128 fmaq(__float128, __float128, __float128);
72 #elif GEOGRAPHICLIB_PRECISION == 5
73 #include <mpreal.h>
74 #endif
75 
76 #if GEOGRAPHICLIB_PRECISION > 3
77 // volatile keyword makes no sense for multiprec types
78 #define GEOGRAPHICLIB_VOLATILE
79 // Signal a convergence failure with multiprec types by throwing an exception
80 // at loop exit.
81 #define GEOGRAPHICLIB_PANIC \
82  (throw GeographicLib::GeographicErr("Convergence failure"), false)
83 #else
84 #define GEOGRAPHICLIB_VOLATILE volatile
85 // Ignore convergence failures with standard floating points types by allowing
86 // loop to exit cleanly.
87 #define GEOGRAPHICLIB_PANIC false
88 #endif
89 
90 namespace GeographicLib {
91 
92  /**
93  * \brief Mathematical functions needed by %GeographicLib
94  *
95  * Define mathematical functions in order to localize system dependencies and
96  * to provide generic versions of the functions. In addition define a real
97  * type to be used by %GeographicLib.
98  *
99  * Example of use:
100  * \include example-Math.cpp
101  **********************************************************************/
103  private:
104  void dummy() {
105  GEOGRAPHICLIB_STATIC_ASSERT(GEOGRAPHICLIB_PRECISION >= 1 &&
107  "Bad value of precision");
108  }
109  Math(); // Disable constructor
110  public:
111 
112 #if GEOGRAPHICLIB_HAVE_LONG_DOUBLE
113  /**
114  * The extended precision type for real numbers, used for some testing.
115  * This is long double on computers with this type; otherwise it is double.
116  **********************************************************************/
117  typedef long double extended;
118 #else
119  typedef double extended;
120 #endif
121 
122 #if GEOGRAPHICLIB_PRECISION == 2
123  /**
124  * The real type for %GeographicLib. Nearly all the testing has been done
125  * with \e real = double. However, the algorithms should also work with
126  * float and long double (where available). (<b>CAUTION</b>: reasonable
127  * accuracy typically cannot be obtained using floats.)
128  **********************************************************************/
129  typedef double real;
130 #elif GEOGRAPHICLIB_PRECISION == 1
131  typedef float real;
132 #elif GEOGRAPHICLIB_PRECISION == 3
133  typedef extended real;
134 #elif GEOGRAPHICLIB_PRECISION == 4
135  typedef boost::multiprecision::float128 real;
136 #elif GEOGRAPHICLIB_PRECISION == 5
137  typedef mpfr::mpreal real;
138 #else
139  typedef double real;
140 #endif
141 
142  /**
143  * @return the number of bits of precision in a real number.
144  **********************************************************************/
145  static inline int digits() {
146 #if GEOGRAPHICLIB_PRECISION != 5
147  return std::numeric_limits<real>::digits;
148 #else
149  return std::numeric_limits<real>::digits();
150 #endif
151  }
152 
153  /**
154  * Set the binary precision of a real number.
155  *
156  * @param[in] ndigits the number of bits of precision.
157  * @return the resulting number of bits of precision.
158  *
159  * This only has an effect when GEOGRAPHICLIB_PRECISION == 5.
160  **********************************************************************/
161  static inline int set_digits(int ndigits) {
162 #if GEOGRAPHICLIB_PRECISION != 5
163  (void)ndigits;
164 #else
165  mpfr::mpreal::set_default_prec(ndigits >= 2 ? ndigits : 2);
166 #endif
167  return digits();
168  }
169 
170  /**
171  * @return the number of decimal digits of precision in a real number.
172  **********************************************************************/
173  static inline int digits10() {
174 #if GEOGRAPHICLIB_PRECISION != 5
175  return std::numeric_limits<real>::digits10;
176 #else
177  return std::numeric_limits<real>::digits10();
178 #endif
179  }
180 
181  /**
182  * Number of additional decimal digits of precision for real relative to
183  * double (0 for float).
184  **********************************************************************/
185  static inline int extra_digits() {
186  return
187  digits10() > std::numeric_limits<double>::digits10 ?
188  digits10() - std::numeric_limits<double>::digits10 : 0;
189  }
190 
191 #if GEOGRAPHICLIB_PRECISION <= 3
192  /**
193  * Number of additional decimal digits of precision of real relative to
194  * double (0 for float).
195  *
196  * <b>DEPRECATED</b>: use extra_digits() instead
197  **********************************************************************/
198  static const int extradigits =
199  std::numeric_limits<real>::digits10 >
200  std::numeric_limits<double>::digits10 ?
201  std::numeric_limits<real>::digits10 -
202  std::numeric_limits<double>::digits10 : 0;
203 #endif
204 
205  /**
206  * true if the machine is big-endian.
207  **********************************************************************/
208  static const bool bigendian = GEOGRAPHICLIB_WORDS_BIGENDIAN;
209 
210  /**
211  * @tparam T the type of the returned value.
212  * @return &pi;.
213  **********************************************************************/
214  template<typename T> static inline T pi() {
215  using std::atan2;
216  static const T pi = atan2(T(0), T(-1));
217  return pi;
218  }
219  /**
220  * A synonym for pi<real>().
221  **********************************************************************/
222  static inline real pi() { return pi<real>(); }
223 
224  /**
225  * @tparam T the type of the returned value.
226  * @return the number of radians in a degree.
227  **********************************************************************/
228  template<typename T> static inline T degree() {
229  static const T degree = pi<T>() / 180;
230  return degree;
231  }
232  /**
233  * A synonym for degree<real>().
234  **********************************************************************/
235  static inline real degree() { return degree<real>(); }
236 
237  /**
238  * Square a number.
239  *
240  * @tparam T the type of the argument and the returned value.
241  * @param[in] x
242  * @return <i>x</i><sup>2</sup>.
243  **********************************************************************/
244  template<typename T> static inline T sq(T x)
245  { return x * x; }
246 
247  /**
248  * The hypotenuse function avoiding underflow and overflow.
249  *
250  * @tparam T the type of the arguments and the returned value.
251  * @param[in] x
252  * @param[in] y
253  * @return sqrt(<i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>).
254  **********************************************************************/
255  template<typename T> static inline T hypot(T x, T y) {
256 #if GEOGRAPHICLIB_CXX11_MATH
257  using std::hypot; return hypot(x, y);
258 #else
259  using std::abs; using std::sqrt;
260  x = abs(x); y = abs(y);
261  if (x < y) std::swap(x, y); // Now x >= y >= 0
262  y /= (x ? x : 1);
263  return x * sqrt(1 + y * y);
264  // For an alternative (square-root free) method see
265  // C. Moler and D. Morrision (1983) https://dx.doi.org/10.1147/rd.276.0577
266  // and A. A. Dubrulle (1983) https://dx.doi.org/10.1147/rd.276.0582
267 #endif
268  }
269 
270  /**
271  * exp(\e x) &minus; 1 accurate near \e x = 0.
272  *
273  * @tparam T the type of the argument and the returned value.
274  * @param[in] x
275  * @return exp(\e x) &minus; 1.
276  **********************************************************************/
277  template<typename T> static inline T expm1(T x) {
278 #if GEOGRAPHICLIB_CXX11_MATH
279  using std::expm1; return expm1(x);
280 #else
281  using std::exp; using std::abs; using std::log;
283  y = exp(x),
284  z = y - 1;
285  // The reasoning here is similar to that for log1p. The expression
286  // mathematically reduces to exp(x) - 1, and the factor z/log(y) = (y -
287  // 1)/log(y) is a slowly varying quantity near y = 1 and is accurately
288  // computed.
289  return abs(x) > 1 ? z : (z == 0 ? x : x * z / log(y));
290 #endif
291  }
292 
293  /**
294  * log(1 + \e x) accurate near \e x = 0.
295  *
296  * @tparam T the type of the argument and the returned value.
297  * @param[in] x
298  * @return log(1 + \e x).
299  **********************************************************************/
300  template<typename T> static inline T log1p(T x) {
301 #if GEOGRAPHICLIB_CXX11_MATH
302  using std::log1p; return log1p(x);
303 #else
304  using std::log;
306  y = 1 + x,
307  z = y - 1;
308  // Here's the explanation for this magic: y = 1 + z, exactly, and z
309  // approx x, thus log(y)/z (which is nearly constant near z = 0) returns
310  // a good approximation to the true log(1 + x)/x. The multiplication x *
311  // (log(y)/z) introduces little additional error.
312  return z == 0 ? x : x * log(y) / z;
313 #endif
314  }
315 
316  /**
317  * The inverse hyperbolic sine function.
318  *
319  * @tparam T the type of the argument and the returned value.
320  * @param[in] x
321  * @return asinh(\e x).
322  **********************************************************************/
323  template<typename T> static inline T asinh(T x) {
324 #if GEOGRAPHICLIB_CXX11_MATH
325  using std::asinh; return asinh(x);
326 #else
327  using std::abs; T y = abs(x); // Enforce odd parity
328  y = log1p(y * (1 + y/(hypot(T(1), y) + 1)));
329  return x < 0 ? -y : y;
330 #endif
331  }
332 
333  /**
334  * The inverse hyperbolic tangent function.
335  *
336  * @tparam T the type of the argument and the returned value.
337  * @param[in] x
338  * @return atanh(\e x).
339  **********************************************************************/
340  template<typename T> static inline T atanh(T x) {
341 #if GEOGRAPHICLIB_CXX11_MATH
342  using std::atanh; return atanh(x);
343 #else
344  using std::abs; T y = abs(x); // Enforce odd parity
345  y = log1p(2 * y/(1 - y))/2;
346  return x < 0 ? -y : y;
347 #endif
348  }
349 
350  /**
351  * The cube root function.
352  *
353  * @tparam T the type of the argument and the returned value.
354  * @param[in] x
355  * @return the real cube root of \e x.
356  **********************************************************************/
357  template<typename T> static inline T cbrt(T x) {
358 #if GEOGRAPHICLIB_CXX11_MATH
359  using std::cbrt; return cbrt(x);
360 #else
361  using std::abs; using std::pow;
362  T y = pow(abs(x), 1/T(3)); // Return the real cube root
363  return x < 0 ? -y : y;
364 #endif
365  }
366 
367  /**
368  * Fused multiply and add.
369  *
370  * @tparam T the type of the arguments and the returned value.
371  * @param[in] x
372  * @param[in] y
373  * @param[in] z
374  * @return <i>xy</i> + <i>z</i>, correctly rounded (on those platforms with
375  * support for the <code>fma</code> instruction).
376  **********************************************************************/
377  template<typename T> static inline T fma(T x, T y, T z) {
378 #if GEOGRAPHICLIB_CXX11_MATH
379  using std::fma; return fma(x, y, z);
380 #else
381  return x * y + z;
382 #endif
383  }
384 
385  /**
386  * Normalize a two-vector.
387  *
388  * @tparam T the type of the argument and the returned value.
389  * @param[in,out] x on output set to <i>x</i>/hypot(<i>x</i>, <i>y</i>).
390  * @param[in,out] y on output set to <i>y</i>/hypot(<i>x</i>, <i>y</i>).
391  **********************************************************************/
392  template<typename T> static inline void norm(T& x, T& y)
393  { T h = hypot(x, y); x /= h; y /= h; }
394 
395  /**
396  * The error-free sum of two numbers.
397  *
398  * @tparam T the type of the argument and the returned value.
399  * @param[in] u
400  * @param[in] v
401  * @param[out] t the exact error given by (\e u + \e v) - \e s.
402  * @return \e s = round(\e u + \e v).
403  *
404  * See D. E. Knuth, TAOCP, Vol 2, 4.2.2, Theorem B. (Note that \e t can be
405  * the same as one of the first two arguments.)
406  **********************************************************************/
407  template<typename T> static inline T sum(T u, T v, T& t) {
408  GEOGRAPHICLIB_VOLATILE T s = u + v;
409  GEOGRAPHICLIB_VOLATILE T up = s - v;
410  GEOGRAPHICLIB_VOLATILE T vpp = s - up;
411  up -= u;
412  vpp -= v;
413  t = -(up + vpp);
414  // u + v = s + t
415  // = round(u + v) + t
416  return s;
417  }
418 
419  /**
420  * Evaluate a polynomial.
421  *
422  * @tparam T the type of the arguments and returned value.
423  * @param[in] N the order of the polynomial.
424  * @param[in] p the coefficient array (of size \e N + 1).
425  * @param[in] x the variable.
426  * @return the value of the polynomial.
427  *
428  * Evaluate <i>y</i> = &sum;<sub><i>n</i>=0..<i>N</i></sub>
429  * <i>p</i><sub><i>n</i></sub> <i>x</i><sup><i>N</i>&minus;<i>n</i></sup>.
430  * Return 0 if \e N &lt; 0. Return <i>p</i><sub>0</sub>, if \e N = 0 (even
431  * if \e x is infinite or a nan). The evaluation uses Horner's method.
432  **********************************************************************/
433  template<typename T> static inline T polyval(int N, const T p[], T x)
434  { T y = N < 0 ? 0 : *p++; while (--N >= 0) y = y * x + *p++; return y; }
435 
436  /**
437  * Normalize an angle (restricted input range).
438  *
439  * @tparam T the type of the argument and returned value.
440  * @param[in] x the angle in degrees.
441  * @return the angle reduced to the range [&minus;180&deg;, 180&deg;).
442  *
443  * \e x must lie in [&minus;540&deg;, 540&deg;).
444  **********************************************************************/
445  template<typename T> static inline T AngNormalize(T x)
446  { return x >= 180 ? x - 360 : (x < -180 ? x + 360 : x); }
447 
448  /**
449  * Normalize an arbitrary angle.
450  *
451  * @tparam T the type of the argument and returned value.
452  * @param[in] x the angle in degrees.
453  * @return the angle reduced to the range [&minus;180&deg;, 180&deg;).
454  *
455  * The range of \e x is unrestricted.
456  **********************************************************************/
457  template<typename T> static inline T AngNormalize2(T x)
458  { using std::fmod; return AngNormalize<T>(fmod(x, T(360))); }
459 
460  /**
461  * Difference of two angles reduced to [&minus;180&deg;, 180&deg;]
462  *
463  * @tparam T the type of the arguments and returned value.
464  * @param[in] x the first angle in degrees.
465  * @param[in] y the second angle in degrees.
466  * @return \e y &minus; \e x, reduced to the range [&minus;180&deg;,
467  * 180&deg;].
468  *
469  * \e x and \e y must both lie in [&minus;180&deg;, 180&deg;]. The result
470  * is equivalent to computing the difference exactly, reducing it to
471  * (&minus;180&deg;, 180&deg;] and rounding the result. Note that this
472  * prescription allows &minus;180&deg; to be returned (e.g., if \e x is
473  * tiny and negative and \e y = 180&deg;).
474  **********************************************************************/
475  template<typename T> static inline T AngDiff(T x, T y) {
476  T t, d = sum(-x, y, t);
477  if ((d - T(180)) + t > T(0)) // y - x > 180
478  d -= T(360); // exact
479  else if ((d + T(180)) + t <= T(0)) // y - x <= -180
480  d += T(360); // exact
481  return d + t;
482  }
483 
484  /**
485  * Coarsen a value close to zero.
486  *
487  * @tparam T the type of the argument and returned value.
488  * @param[in] x
489  * @return the coarsened value.
490  *
491  * The makes the smallest gap in \e x = 1/16 - nextafter(1/16, 0) =
492  * 1/2<sup>57</sup> for reals = 0.7 pm on the earth if \e x is an angle in
493  * degrees. (This is about 1000 times more resolution than we get with
494  * angles around 90&deg;.) We use this to avoid having to deal with near
495  * singular cases when \e x is non-zero but tiny (e.g.,
496  * 10<sup>&minus;200</sup>). This also converts -0 to +0.
497  **********************************************************************/
498  template<typename T> static inline T AngRound(T x) {
499  using std::abs;
500  const T z = 1/T(16);
501  GEOGRAPHICLIB_VOLATILE T y = abs(x);
502  // The compiler mustn't "simplify" z - (z - y) to y
503  y = y < z ? z - (z - y) : y;
504  return x < 0 ? 0 - y : y;
505  }
506 
507  /**
508  * Evaluate the tangent function with the argument in degrees
509  *
510  * @tparam T the type of the argument and the returned value.
511  * @param[in] x in degrees.
512  * @return tan(<i>x</i>).
513  *
514  * If \e x = &plusmn;90&deg;, then a suitably large (but finite) value is
515  * returned.
516  **********************************************************************/
517  template<typename T> static inline T tand(T x) {
518  using std::abs; using std::tan;
519  static const T overflow = 1 / Math::sq(std::numeric_limits<T>::epsilon());
520  return abs(x) != 90 ? tan(x * Math::degree()) :
521  (x < 0 ? -overflow : overflow);
522  }
523 
524  /**
525  * Evaluate the atan function with the result in degrees
526  *
527  * @tparam T the type of the argument and the returned value.
528  * @param[in] x
529  * @return atan(<i>x</i>) in degrees.
530  *
531  * Large values for the argument return &plusmn;90&deg;
532  **********************************************************************/
533  template<typename T> static inline T atand(T x) {
534  using std::abs; using std::atan;
535  static const T
536  overflow = 1 / (Math::sq(std::numeric_limits<T>::epsilon()) * 100);
537  return !(abs(x) >= overflow) ? atan(x) / Math::degree() :
538  (x > 0 ? 90 : -90);
539  }
540 
541  /**
542  * Evaluate the atan2 function with the result in degrees
543  *
544  * @tparam T the type of the arguments and the returned value.
545  * @param[in] y
546  * @param[in] x
547  * @return atan2(<i>y</i>, <i>x</i>) in degrees.
548  *
549  * The result is in the range [&minus;180&deg; 180&deg;).
550  **********************************************************************/
551  template<typename T> static inline T atan2d(T y, T x) {
552  using std::atan2;
553  // The "0 -" converts -0 to +0.
554  return 0 - atan2(-y, x) / Math::degree();
555  }
556 
557  /**
558  * Evaluate <i>e</i> atanh(<i>e x</i>)
559  *
560  * @tparam T the type of the argument and the returned value.
561  * @param[in] x
562  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
563  * sqrt(|<i>e</i><sup>2</sup>|)
564  * @return <i>e</i> atanh(<i>e x</i>)
565  *
566  * If <i>e</i><sup>2</sup> is negative (<i>e</i> is imaginary), the
567  * expression is evaluated in terms of atan.
568  **********************************************************************/
569  template<typename T> static T eatanhe(T x, T es);
570 
571  /**
572  * tan&chi; in terms of tan&phi;
573  *
574  * @tparam T the type of the argument and the returned value.
575  * @param[in] tau &tau; = tan&phi;
576  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
577  * sqrt(|<i>e</i><sup>2</sup>|)
578  * @return &tau;&prime; = tan&chi;
579  *
580  * See Eqs. (7--9) of
581  * C. F. F. Karney,
582  * <a href="https://dx.doi.org/10.1007/s00190-011-0445-3">
583  * Transverse Mercator with an accuracy of a few nanometers,</a>
584  * J. Geodesy 85(8), 475--485 (Aug. 2011)
585  * (preprint <a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>).
586  **********************************************************************/
587  template<typename T> static T taupf(T tau, T es);
588 
589  /**
590  * tan&phi; in terms of tan&chi;
591  *
592  * @tparam T the type of the argument and the returned value.
593  * @param[in] taup &tau;&prime; = tan&chi;
594  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
595  * sqrt(|<i>e</i><sup>2</sup>|)
596  * @return &tau; = tan&phi;
597  *
598  * See Eqs. (19--21) of
599  * C. F. F. Karney,
600  * <a href="https://dx.doi.org/10.1007/s00190-011-0445-3">
601  * Transverse Mercator with an accuracy of a few nanometers,</a>
602  * J. Geodesy 85(8), 475--485 (Aug. 2011)
603  * (preprint <a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>).
604  **********************************************************************/
605  template<typename T> static T tauf(T taup, T es);
606 
607  /**
608  * Test for finiteness.
609  *
610  * @tparam T the type of the argument.
611  * @param[in] x
612  * @return true if number is finite, false if NaN or infinite.
613  **********************************************************************/
614  template<typename T> static inline bool isfinite(T x) {
615 #if GEOGRAPHICLIB_CXX11_MATH
616  using std::isfinite; return isfinite(x);
617 #else
618  using std::abs;
619  return abs(x) <= (std::numeric_limits<T>::max)();
620 #endif
621  }
622 
623  /**
624  * The NaN (not a number)
625  *
626  * @tparam T the type of the returned value.
627  * @return NaN if available, otherwise return the max real of type T.
628  **********************************************************************/
629  template<typename T> static inline T NaN() {
630  return std::numeric_limits<T>::has_quiet_NaN ?
631  std::numeric_limits<T>::quiet_NaN() :
632  (std::numeric_limits<T>::max)();
633  }
634  /**
635  * A synonym for NaN<real>().
636  **********************************************************************/
637  static inline real NaN() { return NaN<real>(); }
638 
639  /**
640  * Test for NaN.
641  *
642  * @tparam T the type of the argument.
643  * @param[in] x
644  * @return true if argument is a NaN.
645  **********************************************************************/
646  template<typename T> static inline bool isnan(T x) {
647 #if GEOGRAPHICLIB_CXX11_MATH
648  using std::isnan; return isnan(x);
649 #else
650  return x != x;
651 #endif
652  }
653 
654  /**
655  * Infinity
656  *
657  * @tparam T the type of the returned value.
658  * @return infinity if available, otherwise return the max real.
659  **********************************************************************/
660  template<typename T> static inline T infinity() {
661  return std::numeric_limits<T>::has_infinity ?
662  std::numeric_limits<T>::infinity() :
663  (std::numeric_limits<T>::max)();
664  }
665  /**
666  * A synonym for infinity<real>().
667  **********************************************************************/
668  static inline real infinity() { return infinity<real>(); }
669 
670  /**
671  * Swap the bytes of a quantity
672  *
673  * @tparam T the type of the argument and the returned value.
674  * @param[in] x
675  * @return x with its bytes swapped.
676  **********************************************************************/
677  template<typename T> static inline T swab(T x) {
678  union {
679  T r;
680  unsigned char c[sizeof(T)];
681  } b;
682  b.r = x;
683  for (int i = sizeof(T)/2; i--; )
684  std::swap(b.c[i], b.c[sizeof(T) - 1 - i]);
685  return b.r;
686  }
687 
688 #if GEOGRAPHICLIB_PRECISION == 4
689  typedef boost::math::policies::policy
690  < boost::math::policies::domain_error
691  <boost::math::policies::errno_on_error>,
692  boost::math::policies::pole_error
693  <boost::math::policies::errno_on_error>,
694  boost::math::policies::overflow_error
695  <boost::math::policies::errno_on_error>,
696  boost::math::policies::evaluation_error
697  <boost::math::policies::errno_on_error> >
698  boost_special_functions_policy;
699 
700  static inline real hypot(real x, real y)
701  { return boost::math::hypot(x, y, boost_special_functions_policy()); }
702 
703  static inline real expm1(real x)
704  { return boost::math::expm1(x, boost_special_functions_policy()); }
705 
706  static inline real log1p(real x)
707  { return boost::math::log1p(x, boost_special_functions_policy()); }
708 
709  static inline real asinh(real x)
710  { return boost::math::asinh(x, boost_special_functions_policy()); }
711 
712  static inline real atanh(real x)
713  { return boost::math::atanh(x, boost_special_functions_policy()); }
714 
715  static inline real cbrt(real x)
716  { return boost::math::cbrt(x, boost_special_functions_policy()); }
717 
718  static inline real fma(real x, real y, real z)
719  { return fmaq(__float128(x), __float128(y), __float128(z)); }
720 
721  static inline bool isnan(real x) { return boost::math::isnan(x); }
722 
723  static inline bool isfinite(real x) { return boost::math::isfinite(x); }
724 #endif
725  };
726 
727 } // namespace GeographicLib
728 
729 #endif // GEOGRAPHICLIB_MATH_HPP
static T AngNormalize(T x)
Definition: Math.hpp:445
static T NaN()
Definition: Math.hpp:629
static T sum(T u, T v, T &t)
Definition: Math.hpp:407
static int digits10()
Definition: Math.hpp:173
static int set_digits(int ndigits)
Definition: Math.hpp:161
static T pi()
Definition: Math.hpp:214
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:90
static T atand(T x)
Definition: Math.hpp:533
static T infinity()
Definition: Math.hpp:660
#define GEOGRAPHICLIB_WORDS_BIGENDIAN
Definition: Math.hpp:40
static real degree()
Definition: Math.hpp:235
GeographicLib::Math::real real
Definition: GeodSolve.cpp:32
static T cbrt(T x)
Definition: Math.hpp:357
static bool isfinite(T x)
Definition: Math.hpp:614
static bool isnan(T x)
Definition: Math.hpp:646
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:102
static T atanh(T x)
Definition: Math.hpp:340
static T expm1(T x)
Definition: Math.hpp:277
static void norm(T &x, T &y)
Definition: Math.hpp:392
#define GEOGRAPHICLIB_PRECISION
Definition: Math.hpp:57
#define GEOGRAPHICLIB_VOLATILE
Definition: Math.hpp:84
static T asinh(T x)
Definition: Math.hpp:323
static int extra_digits()
Definition: Math.hpp:185
static T hypot(T x, T y)
Definition: Math.hpp:255
static T sq(T x)
Definition: Math.hpp:244
static real pi()
Definition: Math.hpp:222
static T atan2d(T y, T x)
Definition: Math.hpp:551
static T polyval(int N, const T p[], T x)
Definition: Math.hpp:433
static T degree()
Definition: Math.hpp:228
static T AngDiff(T x, T y)
Definition: Math.hpp:475
static real NaN()
Definition: Math.hpp:637
static T log1p(T x)
Definition: Math.hpp:300
static T tand(T x)
Definition: Math.hpp:517
static T swab(T x)
Definition: Math.hpp:677
static real infinity()
Definition: Math.hpp:668
Header for GeographicLib::Constants class.
static T fma(T x, T y, T z)
Definition: Math.hpp:377
static T AngRound(T x)
Definition: Math.hpp:498
static int digits()
Definition: Math.hpp:145
static T AngNormalize2(T x)
Definition: Math.hpp:457