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