GeographicLib  1.49
Math.hpp
Go to the documentation of this file.
1 /**
2  * \file Math.hpp
3  * \brief Header for GeographicLib::Math class
4  *
5  * Copyright (c) Charles Karney (2008-2017) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
10 // Constants.hpp includes Math.hpp. Place this include outside Math.hpp's
11 // include guard to enforce this ordering.
13 
14 #if !defined(GEOGRAPHICLIB_MATH_HPP)
15 #define GEOGRAPHICLIB_MATH_HPP 1
16 
17 /**
18  * Are C++11 math functions available?
19  **********************************************************************/
20 #if !defined(GEOGRAPHICLIB_CXX11_MATH)
21 // Recent versions of g++ -std=c++11 (4.7 and later?) set __cplusplus to 201103
22 // and support the new C++11 mathematical functions, std::atanh, etc. However
23 // the Android toolchain, which uses g++ -std=c++11 (4.8 as of 2014-03-11,
24 // according to Pullan Lu), does not support std::atanh. Android toolchains
25 // might define __ANDROID__ or ANDROID; so need to check both. With OSX the
26 // version is GNUC version 4.2 and __cplusplus is set to 201103, so remove the
27 // version check on GNUC.
28 # if defined(__GNUC__) && __cplusplus >= 201103 && \
29  !(defined(__ANDROID__) || defined(ANDROID) || defined(__CYGWIN__))
30 # define GEOGRAPHICLIB_CXX11_MATH 1
31 // Visual C++ 12 supports these functions
32 # elif defined(_MSC_VER) && _MSC_VER >= 1800
33 # define GEOGRAPHICLIB_CXX11_MATH 1
34 # else
35 # define GEOGRAPHICLIB_CXX11_MATH 0
36 # endif
37 #endif
38 
39 #if !defined(GEOGRAPHICLIB_WORDS_BIGENDIAN)
40 # define GEOGRAPHICLIB_WORDS_BIGENDIAN 0
41 #endif
42 
43 #if !defined(GEOGRAPHICLIB_HAVE_LONG_DOUBLE)
44 # define GEOGRAPHICLIB_HAVE_LONG_DOUBLE 0
45 #endif
46 
47 #if !defined(GEOGRAPHICLIB_PRECISION)
48 /**
49  * The precision of floating point numbers used in %GeographicLib. 1 means
50  * float (single precision); 2 (the default) means double; 3 means long double;
51  * 4 is reserved for quadruple precision. Nearly all the testing has been
52  * carried out with doubles and that's the recommended configuration. In order
53  * for long double to be used, GEOGRAPHICLIB_HAVE_LONG_DOUBLE needs to be
54  * defined. Note that with Microsoft Visual Studio, long double is the same as
55  * double.
56  **********************************************************************/
57 # define GEOGRAPHICLIB_PRECISION 2
58 #endif
59 
60 #include <cmath>
61 #include <algorithm>
62 #include <limits>
63 
64 #if GEOGRAPHICLIB_PRECISION == 4
65 #include <boost/version.hpp>
66 #if BOOST_VERSION >= 105600
67 #include <boost/cstdfloat.hpp>
68 #endif
69 #include <boost/multiprecision/float128.hpp>
70 #include <boost/math/special_functions.hpp>
71 __float128 fmaq(__float128, __float128, __float128);
72 #elif GEOGRAPHICLIB_PRECISION == 5
73 #include <mpreal.h>
74 #endif
75 
76 #if GEOGRAPHICLIB_PRECISION > 3
77 // volatile keyword makes no sense for multiprec types
78 #define GEOGRAPHICLIB_VOLATILE
79 // Signal a convergence failure with multiprec types by throwing an exception
80 // at loop exit.
81 #define GEOGRAPHICLIB_PANIC \
82  (throw GeographicLib::GeographicErr("Convergence failure"), false)
83 #else
84 #define GEOGRAPHICLIB_VOLATILE volatile
85 // Ignore convergence failures with standard floating points types by allowing
86 // loop to exit cleanly.
87 #define GEOGRAPHICLIB_PANIC false
88 #endif
89 
90 namespace GeographicLib {
91 
92  /**
93  * \brief Mathematical functions needed by %GeographicLib
94  *
95  * Define mathematical functions in order to localize system dependencies and
96  * to provide generic versions of the functions. In addition define a real
97  * type to be used by %GeographicLib.
98  *
99  * Example of use:
100  * \include example-Math.cpp
101  **********************************************************************/
103  private:
104  void dummy() {
105  GEOGRAPHICLIB_STATIC_ASSERT(GEOGRAPHICLIB_PRECISION >= 1 &&
107  "Bad value of precision");
108  }
109  Math(); // Disable constructor
110  public:
111 
112 #if GEOGRAPHICLIB_HAVE_LONG_DOUBLE
113  /**
114  * The extended precision type for real numbers, used for some testing.
115  * This is long double on computers with this type; otherwise it is double.
116  **********************************************************************/
117  typedef long double extended;
118 #else
119  typedef double extended;
120 #endif
121 
122 #if GEOGRAPHICLIB_PRECISION == 2
123  /**
124  * The real type for %GeographicLib. Nearly all the testing has been done
125  * with \e real = double. However, the algorithms should also work with
126  * float and long double (where available). (<b>CAUTION</b>: reasonable
127  * accuracy typically cannot be obtained using floats.)
128  **********************************************************************/
129  typedef double real;
130 #elif GEOGRAPHICLIB_PRECISION == 1
131  typedef float real;
132 #elif GEOGRAPHICLIB_PRECISION == 3
133  typedef extended real;
134 #elif GEOGRAPHICLIB_PRECISION == 4
135  typedef boost::multiprecision::float128 real;
136 #elif GEOGRAPHICLIB_PRECISION == 5
137  typedef mpfr::mpreal real;
138 #else
139  typedef double real;
140 #endif
141 
142  /**
143  * @return the number of bits of precision in a real number.
144  **********************************************************************/
145  static 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 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 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 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 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 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 T degree() {
217  static const T degree = pi<T>() / 180;
218  return degree;
219  }
220  /**
221  * A synonym for degree<real>().
222  **********************************************************************/
223  static 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 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 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 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 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 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 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 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 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 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 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 T polyval(int N, const T p[], T x)
426  // This used to employ Math::fma; but that's too slow and it seemed not to
427  // improve the accuracy noticeably. This might change when there's direct
428  // hardware support for fma.
429  { T y = N < 0 ? 0 : *p++; while (--N >= 0) y = y * x + *p++; return y; }
430 
431  /**
432  * Normalize an angle.
433  *
434  * @tparam T the type of the argument and returned value.
435  * @param[in] x the angle in degrees.
436  * @return the angle reduced to the range([&minus;180&deg;, 180&deg;].
437  *
438  * The range of \e x is unrestricted.
439  **********************************************************************/
440  template<typename T> static T AngNormalize(T x) {
441 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION != 4
442  using std::remainder;
443  x = remainder(x, T(360)); return x != -180 ? x : 180;
444 #else
445  using std::fmod;
446  T y = fmod(x, T(360));
447 #if defined(_MSC_VER) && _MSC_VER < 1900
448  // Before version 14 (2015), Visual Studio had problems dealing
449  // with -0.0. Specifically
450  // VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0
451  // sincosd has a similar fix.
452  // python 2.7 on Windows 32-bit machines has the same problem.
453  if (x == 0) y = x;
454 #endif
455  return y <= -180 ? y + 360 : (y <= 180 ? y : y - 360);
456 #endif
457  }
458 
459  /**
460  * Normalize a latitude.
461  *
462  * @tparam T the type of the argument and returned value.
463  * @param[in] x the angle in degrees.
464  * @return x if it is in the range [&minus;90&deg;, 90&deg;], otherwise
465  * return NaN.
466  **********************************************************************/
467  template<typename T> static T LatFix(T x)
468  { using std::abs; return abs(x) > 90 ? NaN<T>() : x; }
469 
470  /**
471  * The exact difference of two angles reduced to
472  * (&minus;180&deg;, 180&deg;].
473  *
474  * @tparam T the type of the arguments and returned value.
475  * @param[in] x the first angle in degrees.
476  * @param[in] y the second angle in degrees.
477  * @param[out] e the error term in degrees.
478  * @return \e d, the truncated value of \e y &minus; \e x.
479  *
480  * This computes \e z = \e y &minus; \e x exactly, reduced to
481  * (&minus;180&deg;, 180&deg;]; and then sets \e z = \e d + \e e where \e d
482  * is the nearest representable number to \e z and \e e is the truncation
483  * error. If \e d = &minus;180, then \e e &gt; 0; If \e d = 180, then \e e
484  * &le; 0.
485  **********************************************************************/
486  template<typename T> static T AngDiff(T x, T y, T& e) {
487 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION != 4
488  using std::remainder;
489  T t, d = AngNormalize(sum(remainder(-x, T(360)),
490  remainder( y, T(360)), t));
491 #else
492  T t, d = AngNormalize(sum(AngNormalize(-x), AngNormalize(y), t));
493 #endif
494  // Here y - x = d + t (mod 360), exactly, where d is in (-180,180] and
495  // abs(t) <= eps (eps = 2^-45 for doubles). The only case where the
496  // addition of t takes the result outside the range (-180,180] is d = 180
497  // and t > 0. The case, d = -180 + eps, t = -eps, can't happen, since
498  // sum would have returned the exact result in such a case (i.e., given t
499  // = 0).
500  return sum(d == 180 && t > 0 ? -180 : d, t, e);
501  }
502 
503  /**
504  * Difference of two angles reduced to [&minus;180&deg;, 180&deg;]
505  *
506  * @tparam T the type of the arguments and returned value.
507  * @param[in] x the first angle in degrees.
508  * @param[in] y the second angle in degrees.
509  * @return \e y &minus; \e x, reduced to the range [&minus;180&deg;,
510  * 180&deg;].
511  *
512  * The result is equivalent to computing the difference exactly, reducing
513  * it to (&minus;180&deg;, 180&deg;] and rounding the result. Note that
514  * this prescription allows &minus;180&deg; to be returned (e.g., if \e x
515  * is tiny and negative and \e y = 180&deg;).
516  **********************************************************************/
517  template<typename T> static T AngDiff(T x, T y)
518  { T e; return AngDiff(x, y, e); }
519 
520  /**
521  * Coarsen a value close to zero.
522  *
523  * @tparam T the type of the argument and returned value.
524  * @param[in] x
525  * @return the coarsened value.
526  *
527  * The makes the smallest gap in \e x = 1/16 - nextafter(1/16, 0) =
528  * 1/2<sup>57</sup> for reals = 0.7 pm on the earth if \e x is an angle in
529  * degrees. (This is about 1000 times more resolution than we get with
530  * angles around 90&deg;.) We use this to avoid having to deal with near
531  * singular cases when \e x is non-zero but tiny (e.g.,
532  * 10<sup>&minus;200</sup>). This converts -0 to +0; however tiny negative
533  * numbers get converted to -0.
534  **********************************************************************/
535  template<typename T> static T AngRound(T x) {
536  using std::abs;
537  static const T z = 1/T(16);
538  if (x == 0) return 0;
539  GEOGRAPHICLIB_VOLATILE T y = abs(x);
540  // The compiler mustn't "simplify" z - (z - y) to y
541  y = y < z ? z - (z - y) : y;
542  return x < 0 ? -y : y;
543  }
544 
545  /**
546  * Evaluate the sine and cosine function with the argument in degrees
547  *
548  * @tparam T the type of the arguments.
549  * @param[in] x in degrees.
550  * @param[out] sinx sin(<i>x</i>).
551  * @param[out] cosx cos(<i>x</i>).
552  *
553  * The results obey exactly the elementary properties of the trigonometric
554  * functions, e.g., sin 9&deg; = cos 81&deg; = &minus; sin 123456789&deg;.
555  * If x = &minus;0, then \e sinx = &minus;0; this is the only case where
556  * &minus;0 is returned.
557  **********************************************************************/
558  template<typename T> static void sincosd(T x, T& sinx, T& cosx) {
559  // In order to minimize round-off errors, this function exactly reduces
560  // the argument to the range [-45, 45] before converting it to radians.
561  using std::sin; using std::cos;
562  T r; int q;
563 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION <= 3 && \
564  !defined(__GNUC__)
565  // Disable for gcc because of bug in glibc version < 2.22, see
566  // https://sourceware.org/bugzilla/show_bug.cgi?id=17569
567  // Once this fix is widely deployed, should insert a runtime test for the
568  // glibc version number. For example
569  // #include <gnu/libc-version.h>
570  // std::string version(gnu_get_libc_version()); => "2.22"
571  using std::remquo;
572  r = remquo(x, T(90), &q);
573 #else
574  using std::fmod; using std::floor;
575  r = fmod(x, T(360));
576  q = int(floor(r / 90 + T(0.5)));
577  r -= 90 * q;
578 #endif
579  // now abs(r) <= 45
580  r *= degree();
581  // Possibly could call the gnu extension sincos
582  T s = sin(r), c = cos(r);
583 #if defined(_MSC_VER) && _MSC_VER < 1900
584  // Before version 14 (2015), Visual Studio had problems dealing
585  // with -0.0. Specifically
586  // VC 10,11,12 and 32-bit compile: fmod(-0.0, 360.0) -> +0.0
587  // VC 12 and 64-bit compile: sin(-0.0) -> +0.0
588  // AngNormalize has a similar fix.
589  // python 2.7 on Windows 32-bit machines has the same problem.
590  if (x == 0) s = x;
591 #endif
592  switch (unsigned(q) & 3U) {
593  case 0U: sinx = s; cosx = c; break;
594  case 1U: sinx = c; cosx = -s; break;
595  case 2U: sinx = -s; cosx = -c; break;
596  default: sinx = -c; cosx = s; break; // case 3U
597  }
598  // Set sign of 0 results. -0 only produced for sin(-0)
599  if (x != 0) { sinx += T(0); cosx += T(0); }
600  }
601 
602  /**
603  * Evaluate the sine function with the argument in degrees
604  *
605  * @tparam T the type of the argument and the returned value.
606  * @param[in] x in degrees.
607  * @return sin(<i>x</i>).
608  **********************************************************************/
609  template<typename T> static T sind(T x) {
610  // See sincosd
611  using std::sin; using std::cos;
612  T r; int q;
613 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION <= 3 && \
614  !defined(__GNUC__)
615  using std::remquo;
616  r = remquo(x, T(90), &q);
617 #else
618  using std::fmod; using std::floor;
619  r = fmod(x, T(360));
620  q = int(floor(r / 90 + T(0.5)));
621  r -= 90 * q;
622 #endif
623  // now abs(r) <= 45
624  r *= degree();
625  unsigned p = unsigned(q);
626  r = p & 1U ? cos(r) : sin(r);
627  if (p & 2U) r = -r;
628  if (x != 0) r += T(0);
629  return r;
630  }
631 
632  /**
633  * Evaluate the cosine function with the argument in degrees
634  *
635  * @tparam T the type of the argument and the returned value.
636  * @param[in] x in degrees.
637  * @return cos(<i>x</i>).
638  **********************************************************************/
639  template<typename T> static T cosd(T x) {
640  // See sincosd
641  using std::sin; using std::cos;
642  T r; int q;
643 #if GEOGRAPHICLIB_CXX11_MATH && GEOGRAPHICLIB_PRECISION <= 3 && \
644  !defined(__GNUC__)
645  using std::remquo;
646  r = remquo(x, T(90), &q);
647 #else
648  using std::fmod; using std::floor;
649  r = fmod(x, T(360));
650  q = int(floor(r / 90 + T(0.5)));
651  r -= 90 * q;
652 #endif
653  // now abs(r) <= 45
654  r *= degree();
655  unsigned p = unsigned(q + 1);
656  r = p & 1U ? cos(r) : sin(r);
657  if (p & 2U) r = -r;
658  return T(0) + r;
659  }
660 
661  /**
662  * Evaluate the tangent function with the argument in degrees
663  *
664  * @tparam T the type of the argument and the returned value.
665  * @param[in] x in degrees.
666  * @return tan(<i>x</i>).
667  *
668  * If \e x = &plusmn;90&deg;, then a suitably large (but finite) value is
669  * returned.
670  **********************************************************************/
671  template<typename T> static T tand(T x) {
672  static const T overflow = 1 / sq(std::numeric_limits<T>::epsilon());
673  T s, c;
674  sincosd(x, s, c);
675  return c != 0 ? s / c : (s < 0 ? -overflow : overflow);
676  }
677 
678  /**
679  * Evaluate the atan2 function with the result in degrees
680  *
681  * @tparam T the type of the arguments and the returned value.
682  * @param[in] y
683  * @param[in] x
684  * @return atan2(<i>y</i>, <i>x</i>) in degrees.
685  *
686  * The result is in the range (&minus;180&deg; 180&deg;]. N.B.,
687  * atan2d(&plusmn;0, &minus;1) = +180&deg;; atan2d(&minus;&epsilon;,
688  * &minus;1) = &minus;180&deg;, for &epsilon; positive and tiny;
689  * atan2d(&plusmn;0, +1) = &plusmn;0&deg;.
690  **********************************************************************/
691  template<typename T> static T atan2d(T y, T x) {
692  // In order to minimize round-off errors, this function rearranges the
693  // arguments so that result of atan2 is in the range [-pi/4, pi/4] before
694  // converting it to degrees and mapping the result to the correct
695  // quadrant.
696  using std::atan2; using std::abs;
697  int q = 0;
698  if (abs(y) > abs(x)) { std::swap(x, y); q = 2; }
699  if (x < 0) { x = -x; ++q; }
700  // here x >= 0 and x >= abs(y), so angle is in [-pi/4, pi/4]
701  T ang = atan2(y, x) / degree();
702  switch (q) {
703  // Note that atan2d(-0.0, 1.0) will return -0. However, we expect that
704  // atan2d will not be called with y = -0. If need be, include
705  //
706  // case 0: ang = 0 + ang; break;
707  //
708  // and handle mpfr as in AngRound.
709  case 1: ang = (y >= 0 ? 180 : -180) - ang; break;
710  case 2: ang = 90 - ang; break;
711  case 3: ang = -90 + ang; break;
712  }
713  return ang;
714  }
715 
716  /**
717  * Evaluate the atan function with the result in degrees
718  *
719  * @tparam T the type of the argument and the returned value.
720  * @param[in] x
721  * @return atan(<i>x</i>) in degrees.
722  **********************************************************************/
723  template<typename T> static T atand(T x)
724  { return atan2d(x, T(1)); }
725 
726  /**
727  * Evaluate <i>e</i> atanh(<i>e x</i>)
728  *
729  * @tparam T the type of the argument and the returned value.
730  * @param[in] x
731  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
732  * sqrt(|<i>e</i><sup>2</sup>|)
733  * @return <i>e</i> atanh(<i>e x</i>)
734  *
735  * If <i>e</i><sup>2</sup> is negative (<i>e</i> is imaginary), the
736  * expression is evaluated in terms of atan.
737  **********************************************************************/
738  template<typename T> static T eatanhe(T x, T es);
739 
740  /**
741  * Copy the sign.
742  *
743  * @tparam T the type of the argument.
744  * @param[in] x gives the magitude of the result.
745  * @param[in] y gives the sign of the result.
746  * @return value with the magnitude of \e x and with the sign of \e y.
747  *
748  * This routine correctly handles the case \e y = &minus;0, returning
749  * &minus|<i>x</i>|.
750  **********************************************************************/
751  template<typename T> static T copysign(T x, T y) {
752 #if GEOGRAPHICLIB_CXX11_MATH
753  using std::copysign; return copysign(x, y);
754 #else
755  using std::abs;
756  // NaN counts as positive
757  return abs(x) * (y < 0 || (y == 0 && 1/y < 0) ? -1 : 1);
758 #endif
759  }
760 
761  /**
762  * tan&chi; in terms of tan&phi;
763  *
764  * @tparam T the type of the argument and the returned value.
765  * @param[in] tau &tau; = tan&phi;
766  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
767  * sqrt(|<i>e</i><sup>2</sup>|)
768  * @return &tau;&prime; = tan&chi;
769  *
770  * See Eqs. (7--9) of
771  * C. F. F. Karney,
772  * <a href="https://doi.org/10.1007/s00190-011-0445-3">
773  * Transverse Mercator with an accuracy of a few nanometers,</a>
774  * J. Geodesy 85(8), 475--485 (Aug. 2011)
775  * (preprint
776  * <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>).
777  **********************************************************************/
778  template<typename T> static T taupf(T tau, T es);
779 
780  /**
781  * tan&phi; in terms of tan&chi;
782  *
783  * @tparam T the type of the argument and the returned value.
784  * @param[in] taup &tau;&prime; = tan&chi;
785  * @param[in] es the signed eccentricity = sign(<i>e</i><sup>2</sup>)
786  * sqrt(|<i>e</i><sup>2</sup>|)
787  * @return &tau; = tan&phi;
788  *
789  * See Eqs. (19--21) of
790  * C. F. F. Karney,
791  * <a href="https://doi.org/10.1007/s00190-011-0445-3">
792  * Transverse Mercator with an accuracy of a few nanometers,</a>
793  * J. Geodesy 85(8), 475--485 (Aug. 2011)
794  * (preprint
795  * <a href="https://arxiv.org/abs/1002.1417">arXiv:1002.1417</a>).
796  **********************************************************************/
797  template<typename T> static T tauf(T taup, T es);
798 
799  /**
800  * Test for finiteness.
801  *
802  * @tparam T the type of the argument.
803  * @param[in] x
804  * @return true if number is finite, false if NaN or infinite.
805  **********************************************************************/
806  template<typename T> static bool isfinite(T x) {
807 #if GEOGRAPHICLIB_CXX11_MATH
808  using std::isfinite; return isfinite(x);
809 #else
810  using std::abs;
811 #if defined(_MSC_VER)
812  return abs(x) <= (std::numeric_limits<T>::max)();
813 #else
814  // There's a problem using MPFR C++ 3.6.3 and g++ -std=c++14 (reported on
815  // 2015-05-04) with the parens around std::numeric_limits<T>::max. Of
816  // course, these parens are only needed to deal with Windows stupidly
817  // defining max as a macro. So don't insert the parens on non-Windows
818  // platforms.
819  return abs(x) <= std::numeric_limits<T>::max();
820 #endif
821 #endif
822  }
823 
824  /**
825  * The NaN (not a number)
826  *
827  * @tparam T the type of the returned value.
828  * @return NaN if available, otherwise return the max real of type T.
829  **********************************************************************/
830  template<typename T> static T NaN() {
831 #if defined(_MSC_VER)
832  return std::numeric_limits<T>::has_quiet_NaN ?
833  std::numeric_limits<T>::quiet_NaN() :
834  (std::numeric_limits<T>::max)();
835 #else
836  return std::numeric_limits<T>::has_quiet_NaN ?
837  std::numeric_limits<T>::quiet_NaN() :
838  std::numeric_limits<T>::max();
839 #endif
840  }
841  /**
842  * A synonym for NaN<real>().
843  **********************************************************************/
844  static real NaN() { return NaN<real>(); }
845 
846  /**
847  * Test for NaN.
848  *
849  * @tparam T the type of the argument.
850  * @param[in] x
851  * @return true if argument is a NaN.
852  **********************************************************************/
853  template<typename T> static bool isnan(T x) {
854 #if GEOGRAPHICLIB_CXX11_MATH
855  using std::isnan; return isnan(x);
856 #else
857  return x != x;
858 #endif
859  }
860 
861  /**
862  * Infinity
863  *
864  * @tparam T the type of the returned value.
865  * @return infinity if available, otherwise return the max real.
866  **********************************************************************/
867  template<typename T> static T infinity() {
868 #if defined(_MSC_VER)
869  return std::numeric_limits<T>::has_infinity ?
870  std::numeric_limits<T>::infinity() :
871  (std::numeric_limits<T>::max)();
872 #else
873  return std::numeric_limits<T>::has_infinity ?
874  std::numeric_limits<T>::infinity() :
875  std::numeric_limits<T>::max();
876 #endif
877  }
878  /**
879  * A synonym for infinity<real>().
880  **********************************************************************/
881  static real infinity() { return infinity<real>(); }
882 
883  /**
884  * Swap the bytes of a quantity
885  *
886  * @tparam T the type of the argument and the returned value.
887  * @param[in] x
888  * @return x with its bytes swapped.
889  **********************************************************************/
890  template<typename T> static T swab(T x) {
891  union {
892  T r;
893  unsigned char c[sizeof(T)];
894  } b;
895  b.r = x;
896  for (int i = sizeof(T)/2; i--; )
897  std::swap(b.c[i], b.c[sizeof(T) - 1 - i]);
898  return b.r;
899  }
900 
901 #if GEOGRAPHICLIB_PRECISION == 4
902  typedef boost::math::policies::policy
903  < boost::math::policies::domain_error
904  <boost::math::policies::errno_on_error>,
905  boost::math::policies::pole_error
906  <boost::math::policies::errno_on_error>,
907  boost::math::policies::overflow_error
908  <boost::math::policies::errno_on_error>,
909  boost::math::policies::evaluation_error
910  <boost::math::policies::errno_on_error> >
911  boost_special_functions_policy;
912 
913  static real hypot(real x, real y)
914  { return boost::math::hypot(x, y, boost_special_functions_policy()); }
915 
916  static real expm1(real x)
917  { return boost::math::expm1(x, boost_special_functions_policy()); }
918 
919  static real log1p(real x)
920  { return boost::math::log1p(x, boost_special_functions_policy()); }
921 
922  static real asinh(real x)
923  { return boost::math::asinh(x, boost_special_functions_policy()); }
924 
925  static real atanh(real x)
926  { return boost::math::atanh(x, boost_special_functions_policy()); }
927 
928  static real cbrt(real x)
929  { return boost::math::cbrt(x, boost_special_functions_policy()); }
930 
931  static real fma(real x, real y, real z)
932  { return fmaq(__float128(x), __float128(y), __float128(z)); }
933 
934  static real copysign(real x, real y)
935  { return boost::math::copysign(x, y); }
936 
937  static bool isnan(real x) { return boost::math::isnan(x); }
938 
939  static bool isfinite(real x) { return boost::math::isfinite(x); }
940 #endif
941  };
942 
943 } // namespace GeographicLib
944 
945 #endif // GEOGRAPHICLIB_MATH_HPP
real
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
GeographicLib::Math::AngNormalize
static T AngNormalize(T x)
Definition: Math.hpp:440
GeographicLib::Math::extra_digits
static int extra_digits()
Definition: Math.hpp:187
GeographicLib::Math::sind
static T sind(T x)
Definition: Math.hpp:609
GeographicLib::Math::atanh
static T atanh(T x)
Definition: Math.hpp:328
GeographicLib::Math::NaN
static T NaN()
Definition: Math.hpp:830
GeographicLib::Math::atand
static T atand(T x)
Definition: Math.hpp:723
GeographicLib::Math::tand
static T tand(T x)
Definition: Math.hpp:671
GeographicLib::Math::copysign
static T copysign(T x, T y)
Definition: Math.hpp:751
GeographicLib
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
GeographicLib::Math::infinity
static T infinity()
Definition: Math.hpp:867
GeographicLib::Math::isfinite
static bool isfinite(T x)
Definition: Math.hpp:806
GeographicLib::Math::atan2d
static T atan2d(T y, T x)
Definition: Math.hpp:691
GeographicLib::Math::infinity
static real infinity()
Definition: Math.hpp:881
GeographicLib::Math::sincosd
static void sincosd(T x, T &sinx, T &cosx)
Definition: Math.hpp:558
GeographicLib::Math::extended
double extended
Definition: Math.hpp:119
GeographicLib::Math::expm1
static T expm1(T x)
Definition: Math.hpp:265
GEOGRAPHICLIB_EXPORT
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:91
GeographicLib::Math::pi
static real pi()
Definition: Math.hpp:210
GeographicLib::Math::set_digits
static int set_digits(int ndigits)
Definition: Math.hpp:163
GeographicLib::Math::norm
static void norm(T &x, T &y)
Definition: Math.hpp:384
GeographicLib::Math::polyval
static T polyval(int N, const T p[], T x)
Definition: Math.hpp:425
GeographicLib::Math::real
double real
Definition: Math.hpp:129
GeographicLib::Math
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:102
GeographicLib::Math::degree
static T degree()
Definition: Math.hpp:216
GeographicLib::Math::log1p
static T log1p(T x)
Definition: Math.hpp:288
GeographicLib::Math::AngDiff
static T AngDiff(T x, T y, T &e)
Definition: Math.hpp:486
std::swap
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Definition: NearestNeighbor.hpp:827
GeographicLib::Math::asinh
static T asinh(T x)
Definition: Math.hpp:311
Constants.hpp
Header for GeographicLib::Constants class.
GeographicLib::Math::LatFix
static T LatFix(T x)
Definition: Math.hpp:467
GeographicLib::Math::isnan
static bool isnan(T x)
Definition: Math.hpp:853
GeographicLib::Math::NaN
static real NaN()
Definition: Math.hpp:844
GeographicLib::Math::degree
static real degree()
Definition: Math.hpp:223
GeographicLib::Math::cbrt
static T cbrt(T x)
Definition: Math.hpp:345
GEOGRAPHICLIB_WORDS_BIGENDIAN
#define GEOGRAPHICLIB_WORDS_BIGENDIAN
Definition: Math.hpp:40
GeographicLib::Math::pi
static T pi()
Definition: Math.hpp:202
GeographicLib::Math::cosd
static T cosd(T x)
Definition: Math.hpp:639
GeographicLib::Math::digits
static int digits()
Definition: Math.hpp:145
GEOGRAPHICLIB_PRECISION
#define GEOGRAPHICLIB_PRECISION
Definition: Math.hpp:57
GeographicLib::Math::sum
static T sum(T u, T v, T &t)
Definition: Math.hpp:399
GeographicLib::Math::digits10
static int digits10()
Definition: Math.hpp:175
GeographicLib::Math::AngDiff
static T AngDiff(T x, T y)
Definition: Math.hpp:517
GeographicLib::Math::fma
static T fma(T x, T y, T z)
Definition: Math.hpp:369
GEOGRAPHICLIB_VOLATILE
#define GEOGRAPHICLIB_VOLATILE
Definition: Math.hpp:84
GeographicLib::Math::sq
static T sq(T x)
Definition: Math.hpp:232
GeographicLib::Math::hypot
static T hypot(T x, T y)
Definition: Math.hpp:243
GeographicLib::Math::swab
static T swab(T x)
Definition: Math.hpp:890
GeographicLib::Math::AngRound
static T AngRound(T x)
Definition: Math.hpp:535