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