GeographicLib  1.48
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-2017) <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 /**
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. See also
160  * Utility::set_digits for caveats about when this routine should be
161  * called.
162  **********************************************************************/
163  static inline int set_digits(int ndigits) {
164 #if GEOGRAPHICLIB_PRECISION != 5
165  (void)ndigits;
166 #else
167  mpfr::mpreal::set_default_prec(ndigits >= 2 ? ndigits : 2);
168 #endif
169  return digits();
170  }
171 
172  /**
173  * @return the number of decimal digits of precision in a real number.
174  **********************************************************************/
175  static inline int digits10() {
176 #if GEOGRAPHICLIB_PRECISION != 5
177  return std::numeric_limits<real>::digits10;
178 #else
179  return std::numeric_limits<real>::digits10();
180 #endif
181  }
182 
183  /**
184  * Number of additional decimal digits of precision for real relative to
185  * double (0 for float).
186  **********************************************************************/
187  static inline int extra_digits() {
188  return
189  digits10() > std::numeric_limits<double>::digits10 ?
190  digits10() - std::numeric_limits<double>::digits10 : 0;
191  }
192 
193  /**
194  * true if the machine is big-endian.
195  **********************************************************************/
196  static const bool bigendian = GEOGRAPHICLIB_WORDS_BIGENDIAN;
197 
198  /**
199  * @tparam T the type of the returned value.
200  * @return &pi;.
201  **********************************************************************/
202  template<typename T> static inline T pi() {
203  using std::atan2;
204  static const T pi = atan2(T(0), T(-1));
205  return pi;
206  }
207  /**
208  * A synonym for pi<real>().
209  **********************************************************************/
210  static inline real pi() { return pi<real>(); }
211 
212  /**
213  * @tparam T the type of the returned value.
214  * @return the number of radians in a degree.
215  **********************************************************************/
216  template<typename T> static inline T degree() {
217  static const T degree = pi<T>() / 180;
218  return degree;
219  }
220  /**
221  * A synonym for degree<real>().
222  **********************************************************************/
223  static inline real degree() { return degree<real>(); }
224 
225  /**
226  * Square a number.
227  *
228  * @tparam T the type of the argument and the returned value.
229  * @param[in] x
230  * @return <i>x</i><sup>2</sup>.
231  **********************************************************************/
232  template<typename T> static inline T sq(T x)
233  { return x * x; }
234 
235  /**
236  * The hypotenuse function avoiding underflow and overflow.
237  *
238  * @tparam T the type of the arguments and the returned value.
239  * @param[in] x
240  * @param[in] y
241  * @return sqrt(<i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>).
242  **********************************************************************/
243  template<typename T> static inline T hypot(T x, T y) {
244 #if GEOGRAPHICLIB_CXX11_MATH
245  using std::hypot; return hypot(x, y);
246 #else
247  using std::abs; using std::sqrt;
248  x = abs(x); y = abs(y);
249  if (x < y) std::swap(x, y); // Now x >= y >= 0
250  y /= (x ? x : 1);
251  return x * sqrt(1 + y * y);
252  // For an alternative (square-root free) method see
253  // C. Moler and D. Morrision (1983) https://doi.org/10.1147/rd.276.0577
254  // and A. A. Dubrulle (1983) https://doi.org/10.1147/rd.276.0582
255 #endif
256  }
257 
258  /**
259  * exp(\e x) &minus; 1 accurate near \e x = 0.
260  *
261  * @tparam T the type of the argument and the returned value.
262  * @param[in] x
263  * @return exp(\e x) &minus; 1.
264  **********************************************************************/
265  template<typename T> static inline T expm1(T x) {
266 #if GEOGRAPHICLIB_CXX11_MATH
267  using std::expm1; return expm1(x);
268 #else
269  using std::exp; using std::abs; using std::log;
271  y = exp(x),
272  z = y - 1;
273  // The reasoning here is similar to that for log1p. The expression
274  // mathematically reduces to exp(x) - 1, and the factor z/log(y) = (y -
275  // 1)/log(y) is a slowly varying quantity near y = 1 and is accurately
276  // computed.
277  return abs(x) > 1 ? z : (z == 0 ? x : x * z / log(y));
278 #endif
279  }
280 
281  /**
282  * log(1 + \e x) accurate near \e x = 0.
283  *
284  * @tparam T the type of the argument and the returned value.
285  * @param[in] x
286  * @return log(1 + \e x).
287  **********************************************************************/
288  template<typename T> static inline T log1p(T x) {
289 #if GEOGRAPHICLIB_CXX11_MATH
290  using std::log1p; return log1p(x);
291 #else
292  using std::log;
294  y = 1 + x,
295  z = y - 1;
296  // Here's the explanation for this magic: y = 1 + z, exactly, and z
297  // approx x, thus log(y)/z (which is nearly constant near z = 0) returns
298  // a good approximation to the true log(1 + x)/x. The multiplication x *
299  // (log(y)/z) introduces little additional error.
300  return z == 0 ? x : x * log(y) / z;
301 #endif
302  }
303 
304  /**
305  * The inverse hyperbolic sine function.
306  *
307  * @tparam T the type of the argument and the returned value.
308  * @param[in] x
309  * @return asinh(\e x).
310  **********************************************************************/
311  template<typename T> static inline T asinh(T x) {
312 #if GEOGRAPHICLIB_CXX11_MATH
313  using std::asinh; return asinh(x);
314 #else
315  using std::abs; T y = abs(x); // Enforce odd parity
316  y = log1p(y * (1 + y/(hypot(T(1), y) + 1)));
317  return x < 0 ? -y : y;
318 #endif
319  }
320 
321  /**
322  * The inverse hyperbolic tangent function.
323  *
324  * @tparam T the type of the argument and the returned value.
325  * @param[in] x
326  * @return atanh(\e x).
327  **********************************************************************/
328  template<typename T> static inline T atanh(T x) {
329 #if GEOGRAPHICLIB_CXX11_MATH
330  using std::atanh; return atanh(x);
331 #else
332  using std::abs; T y = abs(x); // Enforce odd parity
333  y = log1p(2 * y/(1 - y))/2;
334  return x < 0 ? -y : y;
335 #endif
336  }
337 
338  /**
339  * The cube root function.
340  *
341  * @tparam T the type of the argument and the returned value.
342  * @param[in] x
343  * @return the real cube root of \e x.
344  **********************************************************************/
345  template<typename T> static inline T cbrt(T x) {
346 #if GEOGRAPHICLIB_CXX11_MATH
347  using std::cbrt; return cbrt(x);
348 #else
349  using std::abs; using std::pow;
350  T y = pow(abs(x), 1/T(3)); // Return the real cube root
351  return x < 0 ? -y : y;
352 #endif
353  }
354 
355  /**
356  * Fused multiply and add.
357  *
358  * @tparam T the type of the arguments and the returned value.
359  * @param[in] x
360  * @param[in] y
361  * @param[in] z
362  * @return <i>xy</i> + <i>z</i>, correctly rounded (on those platforms with
363  * support for the <code>fma</code> instruction).
364  *
365  * On platforms without the <code>fma</code> instruction, no attempt is
366  * made to improve on the result of a rounded multiplication followed by a
367  * rounded addition.
368  **********************************************************************/
369  template<typename T> static inline T fma(T x, T y, T z) {
370 #if GEOGRAPHICLIB_CXX11_MATH
371  using std::fma; return fma(x, y, z);
372 #else
373  return x * y + z;
374 #endif
375  }
376 
377  /**
378  * Normalize a two-vector.
379  *
380  * @tparam T the type of the argument and the returned value.
381  * @param[in,out] x on output set to <i>x</i>/hypot(<i>x</i>, <i>y</i>).
382  * @param[in,out] y on output set to <i>y</i>/hypot(<i>x</i>, <i>y</i>).
383  **********************************************************************/
384  template<typename T> static inline void norm(T& x, T& y)
385  { T h = hypot(x, y); x /= h; y /= h; }
386 
387  /**
388  * The error-free sum of two numbers.
389  *
390  * @tparam T the type of the argument and the returned value.
391  * @param[in] u
392  * @param[in] v
393  * @param[out] t the exact error given by (\e u + \e v) - \e s.
394  * @return \e s = round(\e u + \e v).
395  *
396  * See D. E. Knuth, TAOCP, Vol 2, 4.2.2, Theorem B. (Note that \e t can be
397  * the same as one of the first two arguments.)
398  **********************************************************************/
399  template<typename T> static inline T sum(T u, T v, T& t) {
400  GEOGRAPHICLIB_VOLATILE T s = u + v;
401  GEOGRAPHICLIB_VOLATILE T up = s - v;
402  GEOGRAPHICLIB_VOLATILE T vpp = s - up;
403  up -= u;
404  vpp -= v;
405  t = -(up + vpp);
406  // u + v = s + t
407  // = round(u + v) + t
408  return s;
409  }
410 
411  /**
412  * Evaluate a polynomial.
413  *
414  * @tparam T the type of the arguments and returned value.
415  * @param[in] N the order of the polynomial.
416  * @param[in] p the coefficient array (of size \e N + 1).
417  * @param[in] x the variable.
418  * @return the value of the polynomial.
419  *
420  * Evaluate <i>y</i> = &sum;<sub><i>n</i>=0..<i>N</i></sub>
421  * <i>p</i><sub><i>n</i></sub> <i>x</i><sup><i>N</i>&minus;<i>n</i></sup>.
422  * Return 0 if \e N &lt; 0. Return <i>p</i><sub>0</sub>, if \e N = 0 (even
423  * if \e x is infinite or a nan). The evaluation uses Horner's method.
424  **********************************************************************/
425  template<typename T> static inline T polyval(int N, const T p[], T x)
426  { T y = N < 0 ? 0 : *p++; while (--N >= 0) y = y * x + *p++; return y; }
427 
428  /**
429  * Normalize an angle.
430  *
431  * @tparam T the type of the argument and returned value.
432  * @param[in] x the angle in degrees.
433  * @return the angle reduced to the range([&minus;180&deg;, 180&deg;].
434  *
435  * The range of \e x is unrestricted.
436  **********************************************************************/
437  template<typename T> static inline T AngNormalize(T x) {
438 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION != 4
439  using std::remainder;
440  x = remainder(x, T(360)); return x != -180 ? x : 180;
441 #else
442  using std::fmod;
443  T y = fmod(x, T(360));
444 #if defined(_MSC_VER) && _MSC_VER < 1900
445  // Before version 14 (2015), Visual Studio had problems dealing
446  // with -0.0. Specifically
447  // VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0
448  // sincosd has a similar fix.
449  // python 2.7 on Windows 32-bit machines has the same problem.
450  if (x == 0) y = x;
451 #endif
452  return y <= -180 ? y + 360 : (y <= 180 ? y : y - 360);
453 #endif
454  }
455 
456  /**
457  * Normalize a latitude.
458  *
459  * @tparam T the type of the argument and returned value.
460  * @param[in] x the angle in degrees.
461  * @return x if it is in the range [&minus;90&deg;, 90&deg;], otherwise
462  * return NaN.
463  **********************************************************************/
464  template<typename T> static inline T LatFix(T x)
465  { using std::abs; return abs(x) > 90 ? NaN<T>() : x; }
466 
467  /**
468  * The exact difference of two angles reduced to
469  * (&minus;180&deg;, 180&deg;].
470  *
471  * @tparam T the type of the arguments and returned value.
472  * @param[in] x the first angle in degrees.
473  * @param[in] y the second angle in degrees.
474  * @param[out] e the error term in degrees.
475  * @return \e d, the truncated value of \e y &minus; \e x.
476  *
477  * This computes \e z = \e y &minus; \e x exactly, reduced to
478  * (&minus;180&deg;, 180&deg;]; and then sets \e z = \e d + \e e where \e d
479  * is the nearest representable number to \e z and \e e is the truncation
480  * error. If \e d = &minus;180, then \e e &gt; 0; If \e d = 180, then \e e
481  * &le; 0.
482  **********************************************************************/
483  template<typename T> static inline T AngDiff(T x, T y, T& e) {
484 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION != 4
485  using std::remainder;
486  T t, d = AngNormalize(sum(remainder(-x, T(360)),
487  remainder( y, T(360)), t));
488 #else
489  T t, d = AngNormalize(sum(AngNormalize(-x), AngNormalize(y), t));
490 #endif
491  // Here y - x = d + t (mod 360), exactly, where d is in (-180,180] and
492  // abs(t) <= eps (eps = 2^-45 for doubles). The only case where the
493  // addition of t takes the result outside the range (-180,180] is d = 180
494  // and t > 0. The case, d = -180 + eps, t = -eps, can't happen, since
495  // sum would have returned the exact result in such a case (i.e., given t
496  // = 0).
497  return sum(d == 180 && t > 0 ? -180 : d, t, e);
498  }
499 
500  /**
501  * Difference of two angles reduced to [&minus;180&deg;, 180&deg;]
502  *
503  * @tparam T the type of the arguments and returned value.
504  * @param[in] x the first angle in degrees.
505  * @param[in] y the second angle in degrees.
506  * @return \e y &minus; \e x, reduced to the range [&minus;180&deg;,
507  * 180&deg;].
508  *
509  * The result is equivalent to computing the difference exactly, reducing
510  * it to (&minus;180&deg;, 180&deg;] and rounding the result. Note that
511  * this prescription allows &minus;180&deg; to be returned (e.g., if \e x
512  * is tiny and negative and \e y = 180&deg;).
513  **********************************************************************/
514  template<typename T> static inline T AngDiff(T x, T y)
515  { T e; return AngDiff(x, y, e); }
516 
517  /**
518  * Coarsen a value close to zero.
519  *
520  * @tparam T the type of the argument and returned value.
521  * @param[in] x
522  * @return the coarsened value.
523  *
524  * The makes the smallest gap in \e x = 1/16 - nextafter(1/16, 0) =
525  * 1/2<sup>57</sup> for reals = 0.7 pm on the earth if \e x is an angle in
526  * degrees. (This is about 1000 times more resolution than we get with
527  * angles around 90&deg;.) We use this to avoid having to deal with near
528  * singular cases when \e x is non-zero but tiny (e.g.,
529  * 10<sup>&minus;200</sup>). This converts -0 to +0; however tiny negative
530  * numbers get converted to -0.
531  **********************************************************************/
532  template<typename T> static inline T AngRound(T x) {
533  using std::abs;
534  static const T z = 1/T(16);
535  if (x == 0) return 0;
536  GEOGRAPHICLIB_VOLATILE T y = abs(x);
537  // The compiler mustn't "simplify" z - (z - y) to y
538  y = y < z ? z - (z - y) : y;
539  return x < 0 ? -y : y;
540  }
541 
542  /**
543  * Evaluate the sine and cosine function with the argument in degrees
544  *
545  * @tparam T the type of the arguments.
546  * @param[in] x in degrees.
547  * @param[out] sinx sin(<i>x</i>).
548  * @param[out] cosx cos(<i>x</i>).
549  *
550  * The results obey exactly the elementary properties of the trigonometric
551  * functions, e.g., sin 9&deg; = cos 81&deg; = &minus; sin 123456789&deg;.
552  * If x = &minus;0, then \e sinx = &minus;0; this is the only case where
553  * &minus;0 is returned.
554  **********************************************************************/
555  template<typename T> static inline void sincosd(T x, T& sinx, T& cosx) {
556  // In order to minimize round-off errors, this function exactly reduces
557  // the argument to the range [-45, 45] before converting it to radians.
558  using std::sin; using std::cos;
559  T r; int q;
560 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION <= 3 && \
561  !defined(__GNUC__)
562  // Disable for gcc because of bug in glibc version < 2.22, see
563  // https://sourceware.org/bugzilla/show_bug.cgi?id=17569
564  // Once this fix is widely deployed, should insert a runtime test for the
565  // glibc version number. For example
566  // #include <gnu/libc-version.h>
567  // std::string version(gnu_get_libc_version()); => "2.22"
568  using std::remquo;
569  r = remquo(x, T(90), &q);
570 #else
571  using std::fmod; using std::floor;
572  r = fmod(x, T(360));
573  q = int(floor(r / 90 + T(0.5)));
574  r -= 90 * q;
575 #endif
576  // now abs(r) <= 45
577  r *= degree();
578  // Possibly could call the gnu extension sincos
579  T s = sin(r), c = cos(r);
580 #if defined(_MSC_VER) && _MSC_VER < 1900
581  // Before version 14 (2015), Visual Studio had problems dealing
582  // with -0.0. Specifically
583  // VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0
584  // VC 12 and 64-bit compile: sin(-0.0) -> +0.0
585  // AngNormalize has a similar fix.
586  // python 2.7 on Windows 32-bit machines has the same problem.
587  if (x == 0) s = x;
588 #endif
589  switch (unsigned(q) & 3U) {
590  case 0U: sinx = s; cosx = c; break;
591  case 1U: sinx = c; cosx = -s; break;
592  case 2U: sinx = -s; cosx = -c; break;
593  default: sinx = -c; cosx = s; break; // case 3U
594  }
595  // Set sign of 0 results. -0 only produced for sin(-0)
596  if (x) { sinx += T(0); cosx += T(0); }
597  }
598 
599  /**
600  * Evaluate the sine function with the argument in degrees
601  *
602  * @tparam T the type of the argument and the returned value.
603  * @param[in] x in degrees.
604  * @return sin(<i>x</i>).
605  **********************************************************************/
606  template<typename T> static inline T sind(T x) {
607  // See sincosd
608  using std::sin; using std::cos;
609  T r; int q;
610 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION <= 3 && \
611  !defined(__GNUC__)
612  using std::remquo;
613  r = remquo(x, T(90), &q);
614 #else
615  using std::fmod; using std::floor;
616  r = fmod(x, T(360));
617  q = int(floor(r / 90 + T(0.5)));
618  r -= 90 * q;
619 #endif
620  // now abs(r) <= 45
621  r *= degree();
622  unsigned p = unsigned(q);
623  r = p & 1U ? cos(r) : sin(r);
624  if (p & 2U) r = -r;
625  if (x) r += T(0);
626  return r;
627  }
628 
629  /**
630  * Evaluate the cosine function with the argument in degrees
631  *
632  * @tparam T the type of the argument and the returned value.
633  * @param[in] x in degrees.
634  * @return cos(<i>x</i>).
635  **********************************************************************/
636  template<typename T> static inline T cosd(T x) {
637  // See sincosd
638  using std::sin; using std::cos;
639  T r; int q;
640 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION <= 3 && \
641  !defined(__GNUC__)
642  using std::remquo;
643  r = remquo(x, T(90), &q);
644 #else
645  using std::fmod; using std::floor;
646  r = fmod(x, T(360));
647  q = int(floor(r / 90 + T(0.5)));
648  r -= 90 * q;
649 #endif
650  // now abs(r) <= 45
651  r *= degree();
652  unsigned p = unsigned(q + 1);
653  r = p & 1U ? cos(r) : sin(r);
654  if (p & 2U) r = -r;
655  return T(0) + r;
656  }
657 
658  /**
659  * Evaluate the tangent function with the argument in degrees
660  *
661  * @tparam T the type of the argument and the returned value.
662  * @param[in] x in degrees.
663  * @return tan(<i>x</i>).
664  *
665  * If \e x = &plusmn;90&deg;, then a suitably large (but finite) value is
666  * returned.
667  **********************************************************************/
668  template<typename T> static inline T tand(T x) {
669  static const T overflow = 1 / sq(std::numeric_limits<T>::epsilon());
670  T s, c;
671  sincosd(x, s, c);
672  return c ? s / c : (s < 0 ? -overflow : overflow);
673  }
674 
675  /**
676  * Evaluate the atan2 function with the result in degrees
677  *
678  * @tparam T the type of the arguments and the returned value.
679  * @param[in] y
680  * @param[in] x
681  * @return atan2(<i>y</i>, <i>x</i>) in degrees.
682  *
683  * The result is in the range (&minus;180&deg; 180&deg;]. N.B.,
684  * atan2d(&plusmn;0, &minus;1) = +180&deg;; atan2d(&minus;&epsilon;,
685  * &minus;1) = &minus;180&deg;, for &epsilon; positive and tiny;
686  * atan2d(&plusmn;0, +1) = &plusmn;0&deg;.
687  **********************************************************************/
688  template<typename T> static inline T atan2d(T y, T x) {
689  // In order to minimize round-off errors, this function rearranges the
690  // arguments so that result of atan2 is in the range [-pi/4, pi/4] before
691  // converting it to degrees and mapping the result to the correct
692  // quadrant.
693  using std::atan2; using std::abs;
694  int q = 0;
695  if (abs(y) > abs(x)) { std::swap(x, y); q = 2; }
696  if (x < 0) { x = -x; ++q; }
697  // here x >= 0 and x >= abs(y), so angle is in [-pi/4, pi/4]
698  T ang = atan2(y, x) / degree();
699  switch (q) {
700  // Note that atan2d(-0.0, 1.0) will return -0. However, we expect that
701  // atan2d will not be called with y = -0. If need be, include
702  //
703  // case 0: ang = 0 + ang; break;
704  //
705  // and handle mpfr as in AngRound.
706  case 1: ang = (y >= 0 ? 180 : -180) - ang; break;
707  case 2: ang = 90 - ang; break;
708  case 3: ang = -90 + ang; break;
709  }
710  return ang;
711  }
712 
713  /**
714  * Evaluate the atan function with the result in degrees
715  *
716  * @tparam T the type of the argument and the returned value.
717  * @param[in] x
718  * @return atan(<i>x</i>) in degrees.
719  **********************************************************************/
720  template<typename T> static inline T atand(T x)
721  { return atan2d(x, T(1)); }
722 
723  /**
724  * Evaluate <i>e</i> atanh(<i>e x</i>)
725  *
726  * @tparam T the type of the argument and the returned value.
727  * @param[in] x
728  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
729  * sqrt(|<i>e</i><sup>2</sup>|)
730  * @return <i>e</i> atanh(<i>e x</i>)
731  *
732  * If <i>e</i><sup>2</sup> is negative (<i>e</i> is imaginary), the
733  * expression is evaluated in terms of atan.
734  **********************************************************************/
735  template<typename T> static T eatanhe(T x, T es);
736 
737  /**
738  * Copy the sign.
739  *
740  * @tparam T the type of the argument.
741  * @param[in] x gives the magitude of the result.
742  * @param[in] y gives the sign of the result.
743  * @return value with the magnitude of \e x and with the sign of \e y.
744  *
745  * This routine correctly handles the case \e y = &minus;0, returning
746  * &minus|<i>x</i>|.
747  **********************************************************************/
748  template<typename T> static inline T copysign(T x, T y) {
749 #if GEOGRAPHICLIB_CXX11_MATH
750  using std::copysign; return copysign(x, y);
751 #else
752  using std::abs;
753  // NaN counts as positive
754  return abs(x) * (y < 0 || (y == 0 && 1/y < 0) ? -1 : 1);
755 #endif
756  }
757 
758  /**
759  * tan&chi; in terms of tan&phi;
760  *
761  * @tparam T the type of the argument and the returned value.
762  * @param[in] tau &tau; = tan&phi;
763  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
764  * sqrt(|<i>e</i><sup>2</sup>|)
765  * @return &tau;&prime; = tan&chi;
766  *
767  * See Eqs. (7--9) of
768  * C. F. F. Karney,
769  * <a href="https://doi.org/10.1007/s00190-011-0445-3">
770  * Transverse Mercator with an accuracy of a few nanometers,</a>
771  * J. Geodesy 85(8), 475--485 (Aug. 2011)
772  * (preprint <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>).
773  **********************************************************************/
774  template<typename T> static T taupf(T tau, T es);
775 
776  /**
777  * tan&phi; in terms of tan&chi;
778  *
779  * @tparam T the type of the argument and the returned value.
780  * @param[in] taup &tau;&prime; = tan&chi;
781  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
782  * sqrt(|<i>e</i><sup>2</sup>|)
783  * @return &tau; = tan&phi;
784  *
785  * See Eqs. (19--21) of
786  * C. F. F. Karney,
787  * <a href="https://doi.org/10.1007/s00190-011-0445-3">
788  * Transverse Mercator with an accuracy of a few nanometers,</a>
789  * J. Geodesy 85(8), 475--485 (Aug. 2011)
790  * (preprint <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>).
791  **********************************************************************/
792  template<typename T> static T tauf(T taup, T es);
793 
794  /**
795  * Test for finiteness.
796  *
797  * @tparam T the type of the argument.
798  * @param[in] x
799  * @return true if number is finite, false if NaN or infinite.
800  **********************************************************************/
801  template<typename T> static inline bool isfinite(T x) {
802 #if GEOGRAPHICLIB_CXX11_MATH
803  using std::isfinite; return isfinite(x);
804 #else
805  using std::abs;
806 #if defined(_MSC_VER)
807  return abs(x) <= (std::numeric_limits<T>::max)();
808 #else
809  // There's a problem using MPFR C++ 3.6.3 and g++ -std=c++14 (reported on
810  // 2015-05-04) with the parens around std::numeric_limits<T>::max. Of
811  // course, these parens are only needed to deal with Windows stupidly
812  // defining max as a macro. So don't insert the parens on non-Windows
813  // platforms.
814  return abs(x) <= std::numeric_limits<T>::max();
815 #endif
816 #endif
817  }
818 
819  /**
820  * The NaN (not a number)
821  *
822  * @tparam T the type of the returned value.
823  * @return NaN if available, otherwise return the max real of type T.
824  **********************************************************************/
825  template<typename T> static inline T NaN() {
826 #if defined(_MSC_VER)
827  return std::numeric_limits<T>::has_quiet_NaN ?
828  std::numeric_limits<T>::quiet_NaN() :
829  (std::numeric_limits<T>::max)();
830 #else
831  return std::numeric_limits<T>::has_quiet_NaN ?
832  std::numeric_limits<T>::quiet_NaN() :
833  std::numeric_limits<T>::max();
834 #endif
835  }
836  /**
837  * A synonym for NaN<real>().
838  **********************************************************************/
839  static inline real NaN() { return NaN<real>(); }
840 
841  /**
842  * Test for NaN.
843  *
844  * @tparam T the type of the argument.
845  * @param[in] x
846  * @return true if argument is a NaN.
847  **********************************************************************/
848  template<typename T> static inline bool isnan(T x) {
849 #if GEOGRAPHICLIB_CXX11_MATH
850  using std::isnan; return isnan(x);
851 #else
852  return x != x;
853 #endif
854  }
855 
856  /**
857  * Infinity
858  *
859  * @tparam T the type of the returned value.
860  * @return infinity if available, otherwise return the max real.
861  **********************************************************************/
862  template<typename T> static inline T infinity() {
863 #if defined(_MSC_VER)
864  return std::numeric_limits<T>::has_infinity ?
865  std::numeric_limits<T>::infinity() :
866  (std::numeric_limits<T>::max)();
867 #else
868  return std::numeric_limits<T>::has_infinity ?
869  std::numeric_limits<T>::infinity() :
870  std::numeric_limits<T>::max();
871 #endif
872  }
873  /**
874  * A synonym for infinity<real>().
875  **********************************************************************/
876  static inline real infinity() { return infinity<real>(); }
877 
878  /**
879  * Swap the bytes of a quantity
880  *
881  * @tparam T the type of the argument and the returned value.
882  * @param[in] x
883  * @return x with its bytes swapped.
884  **********************************************************************/
885  template<typename T> static inline T swab(T x) {
886  union {
887  T r;
888  unsigned char c[sizeof(T)];
889  } b;
890  b.r = x;
891  for (int i = sizeof(T)/2; i--; )
892  std::swap(b.c[i], b.c[sizeof(T) - 1 - i]);
893  return b.r;
894  }
895 
896 #if GEOGRAPHICLIB_PRECISION == 4
897  typedef boost::math::policies::policy
898  < boost::math::policies::domain_error
899  <boost::math::policies::errno_on_error>,
900  boost::math::policies::pole_error
901  <boost::math::policies::errno_on_error>,
902  boost::math::policies::overflow_error
903  <boost::math::policies::errno_on_error>,
904  boost::math::policies::evaluation_error
905  <boost::math::policies::errno_on_error> >
906  boost_special_functions_policy;
907 
908  static inline real hypot(real x, real y)
909  { return boost::math::hypot(x, y, boost_special_functions_policy()); }
910 
911  static inline real expm1(real x)
912  { return boost::math::expm1(x, boost_special_functions_policy()); }
913 
914  static inline real log1p(real x)
915  { return boost::math::log1p(x, boost_special_functions_policy()); }
916 
917  static inline real asinh(real x)
918  { return boost::math::asinh(x, boost_special_functions_policy()); }
919 
920  static inline real atanh(real x)
921  { return boost::math::atanh(x, boost_special_functions_policy()); }
922 
923  static inline real cbrt(real x)
924  { return boost::math::cbrt(x, boost_special_functions_policy()); }
925 
926  static inline real fma(real x, real y, real z)
927  { return fmaq(__float128(x), __float128(y), __float128(z)); }
928 
929  static inline real copysign(real x, real y)
930  { return boost::math::copysign(x, y); }
931 
932  static inline bool isnan(real x) { return boost::math::isnan(x); }
933 
934  static inline bool isfinite(real x) { return boost::math::isfinite(x); }
935 #endif
936  };
937 
938 } // namespace GeographicLib
939 
940 #endif // GEOGRAPHICLIB_MATH_HPP
static T AngNormalize(T x)
Definition: Math.hpp:437
static T NaN()
Definition: Math.hpp:825
static T sum(T u, T v, T &t)
Definition: Math.hpp:399
static int digits10()
Definition: Math.hpp:175
static int set_digits(int ndigits)
Definition: Math.hpp:163
static T pi()
Definition: Math.hpp:202
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:91
static T atand(T x)
Definition: Math.hpp:720
static T infinity()
Definition: Math.hpp:862
#define GEOGRAPHICLIB_WORDS_BIGENDIAN
Definition: Math.hpp:40
static real degree()
Definition: Math.hpp:223
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
static T cbrt(T x)
Definition: Math.hpp:345
static bool isfinite(T x)
Definition: Math.hpp:801
static bool isnan(T x)
Definition: Math.hpp:848
static T LatFix(T x)
Definition: Math.hpp:464
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:102
static T sind(T x)
Definition: Math.hpp:606
static void sincosd(T x, T &sinx, T &cosx)
Definition: Math.hpp:555
static T AngDiff(T x, T y, T &e)
Definition: Math.hpp:483
static T atanh(T x)
Definition: Math.hpp:328
static T expm1(T x)
Definition: Math.hpp:265
static void norm(T &x, T &y)
Definition: Math.hpp:384
#define GEOGRAPHICLIB_PRECISION
Definition: Math.hpp:57
#define GEOGRAPHICLIB_VOLATILE
Definition: Math.hpp:84
static T asinh(T x)
Definition: Math.hpp:311
static int extra_digits()
Definition: Math.hpp:187
static T hypot(T x, T y)
Definition: Math.hpp:243
static T sq(T x)
Definition: Math.hpp:232
static real pi()
Definition: Math.hpp:210
static T cosd(T x)
Definition: Math.hpp:636
static T atan2d(T y, T x)
Definition: Math.hpp:688
static T polyval(int N, const T p[], T x)
Definition: Math.hpp:425
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
static T degree()
Definition: Math.hpp:216
static T AngDiff(T x, T y)
Definition: Math.hpp:514
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
static real NaN()
Definition: Math.hpp:839
static T log1p(T x)
Definition: Math.hpp:288
static T tand(T x)
Definition: Math.hpp:668
static T copysign(T x, T y)
Definition: Math.hpp:748
static T swab(T x)
Definition: Math.hpp:885
static real infinity()
Definition: Math.hpp:876
Header for GeographicLib::Constants class.
static T fma(T x, T y, T z)
Definition: Math.hpp:369
static T AngRound(T x)
Definition: Math.hpp:532
static int digits()
Definition: Math.hpp:145