GeographicLib  1.46
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-2016) <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  /**
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://dx.doi.org/10.1147/rd.276.0577
254  // and A. A. Dubrulle (1983) https://dx.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  x = fmod(x, T(360));
444  return x < -180 ? x + 360 : (x < 180 ? x : x - 360);
445 #endif
446  }
447 
448  /**
449  * Normalize a latitude.
450  *
451  * @tparam T the type of the argument and returned value.
452  * @param[in] x the angle in degrees.
453  * @return x if it is in the range [&minus;90&deg;, 90&deg;], otherwise
454  * return NaN.
455  **********************************************************************/
456  template<typename T> static inline T LatFix(T x)
457  { using std::abs; return abs(x) > 90 ? NaN<T>() : x; }
458 
459  /**
460  * The exact difference of two angles reduced to
461  * (&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  * @param[out] e the error term in degrees.
467  * @return \e d, the truncated value of \e y &minus; \e x.
468  *
469  * This computes \e z = \e y &minus; \e x exactly, reduced to
470  * (&minus;180&deg;, 180&deg;]; and then sets \e z = \e d + \e e where \e d
471  * is the nearest representable number to \e z and \e e is the truncation
472  * error. If \e d = &minus;180, then \e e &gt; 0; If \e d = 180, then \e e
473  * &le; 0.
474  **********************************************************************/
475  template<typename T> static inline T AngDiff(T x, T y, T& e) {
476 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION != 4
477  using std::remainder;
478  T t, d = - AngNormalize(sum(remainder( x, T(360)),
479  remainder(-y, T(360)), t));
480 #else
481  T t, d = - AngNormalize(sum(AngNormalize(x), AngNormalize(-y), t));
482 #endif
483  // Here y - x = d - t (mod 360), exactly, where d is in (-180,180] and
484  // abs(t) <= eps (eps = 2^-45 for doubles). The only case where the
485  // addition of t takes the result outside the range (-180,180] is d = 180
486  // and t < 0. The case, d = -180 + eps, t = eps, can't happen, since
487  // sum would have returned the exact result in such a case (i.e., given t
488  // = 0).
489  return sum(d == 180 && t < 0 ? -180 : d, -t, e);
490  }
491 
492  /**
493  * Difference of two angles reduced to [&minus;180&deg;, 180&deg;]
494  *
495  * @tparam T the type of the arguments and returned value.
496  * @param[in] x the first angle in degrees.
497  * @param[in] y the second angle in degrees.
498  * @return \e y &minus; \e x, reduced to the range [&minus;180&deg;,
499  * 180&deg;].
500  *
501  * The result is equivalent to computing the difference exactly, reducing
502  * it to (&minus;180&deg;, 180&deg;] and rounding the result. Note that
503  * this prescription allows &minus;180&deg; to be returned (e.g., if \e x
504  * is tiny and negative and \e y = 180&deg;).
505  **********************************************************************/
506  template<typename T> static inline T AngDiff(T x, T y)
507  { T e; return AngDiff(x, y, e); }
508 
509  /**
510  * Coarsen a value close to zero.
511  *
512  * @tparam T the type of the argument and returned value.
513  * @param[in] x
514  * @return the coarsened value.
515  *
516  * The makes the smallest gap in \e x = 1/16 - nextafter(1/16, 0) =
517  * 1/2<sup>57</sup> for reals = 0.7 pm on the earth if \e x is an angle in
518  * degrees. (This is about 1000 times more resolution than we get with
519  * angles around 90&deg;.) We use this to avoid having to deal with near
520  * singular cases when \e x is non-zero but tiny (e.g.,
521  * 10<sup>&minus;200</sup>). This converts -0 to +0; however tiny negative
522  * numbers get converted to -0.
523  **********************************************************************/
524  template<typename T> static inline T AngRound(T x) {
525  using std::abs;
526  static const T z = 1/T(16);
527  if (x == 0) return 0;
528  GEOGRAPHICLIB_VOLATILE T y = abs(x);
529  // The compiler mustn't "simplify" z - (z - y) to y
530  y = y < z ? z - (z - y) : y;
531  return x < 0 ? -y : y;
532  }
533 
534  /**
535  * Evaluate the sine and cosine function with the argument in degrees
536  *
537  * @tparam T the type of the arguments.
538  * @param[in] x in degrees.
539  * @param[out] sinx sin(<i>x</i>).
540  * @param[out] cosx cos(<i>x</i>).
541  *
542  * The results obey exactly the elementary properties of the trigonometric
543  * functions, e.g., sin 9&deg; = cos 81&deg; = &minus; sin 123456789&deg;.
544  **********************************************************************/
545  template<typename T> static inline void sincosd(T x, T& sinx, T& cosx) {
546  // In order to minimize round-off errors, this function exactly reduces
547  // the argument to the range [-45, 45] before converting it to radians.
548  using std::sin; using std::cos;
549  T r; int q;
550 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION <= 3 && \
551  !defined(__GNUC__)
552  // Disable for gcc because of bug in glibc version < 2.22, see
553  // https://sourceware.org/bugzilla/show_bug.cgi?id=17569
554  // Once this fix is widely deployed, should insert a runtime test for the
555  // glibc version number. For example
556  // #include <gnu/libc-version.h>
557  // std::string version(gnu_get_libc_version()); => "2.22"
558  using std::remquo;
559  r = remquo(x, T(90), &q);
560 #else
561  using std::fmod; using std::floor;
562  r = fmod(x, T(360));
563  q = int(floor(r / 90 + T(0.5)));
564  r -= 90 * q;
565 #endif
566  // now abs(r) <= 45
567  r *= degree();
568  // Possibly could call the gnu extension sincos
569  T s = sin(r), c = cos(r);
570 #if defined(_MSC_VER) && _MSC_VER < 1900
571  // Before version 14 (2015), Visual Studio had problems dealing
572  // with -0.0. Specifically
573  // VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0
574  // VC 12 and 64-bit compile: sin(-0.0) -> +0.0
575  if (x == 0) s = x;
576 #endif
577  switch (unsigned(q) & 3U) {
578  case 0U: sinx = s; cosx = c; break;
579  case 1U: sinx = c; cosx = T(0) - s; break;
580  case 2U: sinx = T(0) - s; cosx = T(0) - c; break;
581  default: sinx = T(0) - c; cosx = s; break; // case 3U
582  }
583  }
584 
585  /**
586  * Evaluate the sine function with the argument in degrees
587  *
588  * @tparam T the type of the argument and the returned value.
589  * @param[in] x in degrees.
590  * @return sin(<i>x</i>).
591  **********************************************************************/
592  template<typename T> static inline T sind(T x) {
593  // See sincosd
594  using std::sin; using std::cos;
595  T r; int q;
596 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION <= 3 && \
597  !defined(__GNUC__)
598  using std::remquo;
599  r = remquo(x, T(90), &q);
600 #else
601  using std::fmod; using std::floor;
602  r = fmod(x, T(360));
603  q = int(floor(r / 90 + T(0.5)));
604  r -= 90 * q;
605 #endif
606  // now abs(r) <= 45
607  r *= degree();
608  unsigned p = unsigned(q);
609  r = p & 1U ? cos(r) : sin(r);
610  return p & 2U ? T(0) - r : r;
611  }
612 
613  /**
614  * Evaluate the cosine function with the argument in degrees
615  *
616  * @tparam T the type of the argument and the returned value.
617  * @param[in] x in degrees.
618  * @return cos(<i>x</i>).
619  **********************************************************************/
620  template<typename T> static inline T cosd(T x) {
621  // See sincosd
622  using std::sin; using std::cos;
623  T r; int q;
624 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION <= 3 && \
625  !defined(__GNUC__)
626  using std::remquo;
627  r = remquo(x, T(90), &q);
628 #else
629  using std::fmod; using std::floor;
630  r = fmod(x, T(360));
631  q = int(floor(r / 90 + T(0.5)));
632  r -= 90 * q;
633 #endif
634  // now abs(r) <= 45
635  r *= degree();
636  unsigned p = unsigned(q + 1);
637  r = p & 1U ? cos(r) : sin(r);
638  return p & 2U ? T(0) - r : r;
639  }
640 
641  /**
642  * Evaluate the tangent function with the argument in degrees
643  *
644  * @tparam T the type of the argument and the returned value.
645  * @param[in] x in degrees.
646  * @return tan(<i>x</i>).
647  *
648  * If \e x = &plusmn;90&deg;, then a suitably large (but finite) value is
649  * returned.
650  **********************************************************************/
651  template<typename T> static inline T tand(T x) {
652  static const T overflow = 1 / sq(std::numeric_limits<T>::epsilon());
653  T s, c;
654  sincosd(x, s, c);
655  return c ? s / c : (s < 0 ? -overflow : overflow);
656  }
657 
658  /**
659  * Evaluate the atan2 function with the result in degrees
660  *
661  * @tparam T the type of the arguments and the returned value.
662  * @param[in] y
663  * @param[in] x
664  * @return atan2(<i>y</i>, <i>x</i>) in degrees.
665  *
666  * The result is in the range [&minus;180&deg; 180&deg;). N.B.,
667  * atan2d(&plusmn;0, &minus;1) = &minus;180&deg;; atan2d(+&epsilon;,
668  * &minus;1) = +180&deg;, for &epsilon; positive and tiny;
669  * atan2d(&plusmn;0, 1) = &plusmn;0&deg;.
670  **********************************************************************/
671  template<typename T> static inline T atan2d(T y, T x) {
672  // In order to minimize round-off errors, this function rearranges the
673  // arguments so that result of atan2 is in the range [-pi/4, pi/4] before
674  // converting it to degrees and mapping the result to the correct
675  // quadrant.
676  using std::atan2; using std::abs;
677  int q = 0;
678  if (abs(y) > abs(x)) { std::swap(x, y); q = 2; }
679  if (x < 0) { x = -x; ++q; }
680  // here x >= 0 and x >= abs(y), so angle is in [-pi/4, pi/4]
681  T ang = atan2(y, x) / degree();
682  switch (q) {
683  // Note that atan2d(-0.0, 1.0) will return -0. However, we expect that
684  // atan2d will not be called with y = -0. If need be, include
685  //
686  // case 0: ang = 0 + ang; break;
687  //
688  // and handle mpfr as in AngRound.
689  case 1: ang = (y > 0 ? 180 : -180) - ang; break;
690  case 2: ang = 90 - ang; break;
691  case 3: ang = -90 + ang; break;
692  }
693  return ang;
694  }
695 
696  /**
697  * Evaluate the atan function with the result in degrees
698  *
699  * @tparam T the type of the argument and the returned value.
700  * @param[in] x
701  * @return atan(<i>x</i>) in degrees.
702  **********************************************************************/
703  template<typename T> static inline T atand(T x)
704  { return atan2d(x, T(1)); }
705 
706  /**
707  * Evaluate <i>e</i> atanh(<i>e x</i>)
708  *
709  * @tparam T the type of the argument and the returned value.
710  * @param[in] x
711  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
712  * sqrt(|<i>e</i><sup>2</sup>|)
713  * @return <i>e</i> atanh(<i>e x</i>)
714  *
715  * If <i>e</i><sup>2</sup> is negative (<i>e</i> is imaginary), the
716  * expression is evaluated in terms of atan.
717  **********************************************************************/
718  template<typename T> static T eatanhe(T x, T es);
719 
720  /**
721  * Copy the sign.
722  *
723  * @tparam T the type of the argument.
724  * @param[in] x gives the magitude of the result.
725  * @param[in] y gives the sign of the result.
726  * @return value with the magnitude of \e x and with the sign of \e y.
727  **********************************************************************/
728  template<typename T> static inline T copysign(T x, T y) {
729 #if GEOGRAPHICLIB_CXX11_MATH
730  using std::copysign; return copysign(x, y);
731 #else
732  using std::abs; using std::atan2;
733  // NaN counts as positive
734  return abs(x) * (y < 0 || (y == 0 && 1/y < 0) ? -1 : 1);
735 #endif
736  }
737 
738  /**
739  * tan&chi; in terms of tan&phi;
740  *
741  * @tparam T the type of the argument and the returned value.
742  * @param[in] tau &tau; = tan&phi;
743  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
744  * sqrt(|<i>e</i><sup>2</sup>|)
745  * @return &tau;&prime; = tan&chi;
746  *
747  * See Eqs. (7--9) of
748  * C. F. F. Karney,
749  * <a href="https://dx.doi.org/10.1007/s00190-011-0445-3">
750  * Transverse Mercator with an accuracy of a few nanometers,</a>
751  * J. Geodesy 85(8), 475--485 (Aug. 2011)
752  * (preprint <a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>).
753  **********************************************************************/
754  template<typename T> static T taupf(T tau, T es);
755 
756  /**
757  * tan&phi; in terms of tan&chi;
758  *
759  * @tparam T the type of the argument and the returned value.
760  * @param[in] taup &tau;&prime; = tan&chi;
761  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
762  * sqrt(|<i>e</i><sup>2</sup>|)
763  * @return &tau; = tan&phi;
764  *
765  * See Eqs. (19--21) of
766  * C. F. F. Karney,
767  * <a href="https://dx.doi.org/10.1007/s00190-011-0445-3">
768  * Transverse Mercator with an accuracy of a few nanometers,</a>
769  * J. Geodesy 85(8), 475--485 (Aug. 2011)
770  * (preprint <a href="http://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>).
771  **********************************************************************/
772  template<typename T> static T tauf(T taup, T es);
773 
774  /**
775  * Test for finiteness.
776  *
777  * @tparam T the type of the argument.
778  * @param[in] x
779  * @return true if number is finite, false if NaN or infinite.
780  **********************************************************************/
781  template<typename T> static inline bool isfinite(T x) {
782 #if GEOGRAPHICLIB_CXX11_MATH
783  using std::isfinite; return isfinite(x);
784 #else
785  using std::abs;
786 #if defined(_MSC_VER)
787  return abs(x) <= (std::numeric_limits<T>::max)();
788 #else
789  // There's a problem using MPFR C++ 3.6.3 and g++ -std=c++14 (reported on
790  // 2015-05-04) with the parens around std::numeric_limits<T>::max. Of
791  // course, these parens are only needed to deal with Windows stupidly
792  // defining max as a macro. So don't insert the parens on non-Windows
793  // platforms.
794  return abs(x) <= std::numeric_limits<T>::max();
795 #endif
796 #endif
797  }
798 
799  /**
800  * The NaN (not a number)
801  *
802  * @tparam T the type of the returned value.
803  * @return NaN if available, otherwise return the max real of type T.
804  **********************************************************************/
805  template<typename T> static inline T NaN() {
806 #if defined(_MSC_VER)
807  return std::numeric_limits<T>::has_quiet_NaN ?
808  std::numeric_limits<T>::quiet_NaN() :
809  (std::numeric_limits<T>::max)();
810 #else
811  return std::numeric_limits<T>::has_quiet_NaN ?
812  std::numeric_limits<T>::quiet_NaN() :
813  std::numeric_limits<T>::max();
814 #endif
815  }
816  /**
817  * A synonym for NaN<real>().
818  **********************************************************************/
819  static inline real NaN() { return NaN<real>(); }
820 
821  /**
822  * Test for NaN.
823  *
824  * @tparam T the type of the argument.
825  * @param[in] x
826  * @return true if argument is a NaN.
827  **********************************************************************/
828  template<typename T> static inline bool isnan(T x) {
829 #if GEOGRAPHICLIB_CXX11_MATH
830  using std::isnan; return isnan(x);
831 #else
832  return x != x;
833 #endif
834  }
835 
836  /**
837  * Infinity
838  *
839  * @tparam T the type of the returned value.
840  * @return infinity if available, otherwise return the max real.
841  **********************************************************************/
842  template<typename T> static inline T infinity() {
843 #if defined(_MSC_VER)
844  return std::numeric_limits<T>::has_infinity ?
845  std::numeric_limits<T>::infinity() :
846  (std::numeric_limits<T>::max)();
847 #else
848  return std::numeric_limits<T>::has_infinity ?
849  std::numeric_limits<T>::infinity() :
850  std::numeric_limits<T>::max();
851 #endif
852  }
853  /**
854  * A synonym for infinity<real>().
855  **********************************************************************/
856  static inline real infinity() { return infinity<real>(); }
857 
858  /**
859  * Swap the bytes of a quantity
860  *
861  * @tparam T the type of the argument and the returned value.
862  * @param[in] x
863  * @return x with its bytes swapped.
864  **********************************************************************/
865  template<typename T> static inline T swab(T x) {
866  union {
867  T r;
868  unsigned char c[sizeof(T)];
869  } b;
870  b.r = x;
871  for (int i = sizeof(T)/2; i--; )
872  std::swap(b.c[i], b.c[sizeof(T) - 1 - i]);
873  return b.r;
874  }
875 
876 #if GEOGRAPHICLIB_PRECISION == 4
877  typedef boost::math::policies::policy
878  < boost::math::policies::domain_error
879  <boost::math::policies::errno_on_error>,
880  boost::math::policies::pole_error
881  <boost::math::policies::errno_on_error>,
882  boost::math::policies::overflow_error
883  <boost::math::policies::errno_on_error>,
884  boost::math::policies::evaluation_error
885  <boost::math::policies::errno_on_error> >
886  boost_special_functions_policy;
887 
888  static inline real hypot(real x, real y)
889  { return boost::math::hypot(x, y, boost_special_functions_policy()); }
890 
891  static inline real expm1(real x)
892  { return boost::math::expm1(x, boost_special_functions_policy()); }
893 
894  static inline real log1p(real x)
895  { return boost::math::log1p(x, boost_special_functions_policy()); }
896 
897  static inline real asinh(real x)
898  { return boost::math::asinh(x, boost_special_functions_policy()); }
899 
900  static inline real atanh(real x)
901  { return boost::math::atanh(x, boost_special_functions_policy()); }
902 
903  static inline real cbrt(real x)
904  { return boost::math::cbrt(x, boost_special_functions_policy()); }
905 
906  static inline real fma(real x, real y, real z)
907  { return fmaq(__float128(x), __float128(y), __float128(z)); }
908 
909  static inline real copysign(real x, real y)
910  { return boost::math::copysign(x, y); }
911 
912  static inline bool isnan(real x) { return boost::math::isnan(x); }
913 
914  static inline bool isfinite(real x) { return boost::math::isfinite(x); }
915 #endif
916  };
917 
918 } // namespace GeographicLib
919 
920 #endif // GEOGRAPHICLIB_MATH_HPP
static T AngNormalize(T x)
Definition: Math.hpp:437
static T NaN()
Definition: Math.hpp:805
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:90
static T atand(T x)
Definition: Math.hpp:703
static T infinity()
Definition: Math.hpp:842
#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:781
static bool isnan(T x)
Definition: Math.hpp:828
static T LatFix(T x)
Definition: Math.hpp:456
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:102
static T sind(T x)
Definition: Math.hpp:592
static void sincosd(T x, T &sinx, T &cosx)
Definition: Math.hpp:545
static T AngDiff(T x, T y, T &e)
Definition: Math.hpp:475
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:620
static T atan2d(T y, T x)
Definition: Math.hpp:671
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:506
static real NaN()
Definition: Math.hpp:819
static T log1p(T x)
Definition: Math.hpp:288
static T tand(T x)
Definition: Math.hpp:651
static T copysign(T x, T y)
Definition: Math.hpp:728
static T swab(T x)
Definition: Math.hpp:865
static real infinity()
Definition: Math.hpp:856
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:524
static int digits()
Definition: Math.hpp:145