GeographicLib  1.38
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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-2014) <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.
26 # if defined(__GNUC__) && __GNUC__ == 4 && __GNUC_MINOR__ >= 7 \
27  && __cplusplus >= 201103 && \
28  !(defined(__ANDROID__) || defined(ANDROID) || defined(__CYGWIN__))
29 # define GEOGRAPHICLIB_CXX11_MATH 1
30 // Visual C++ 12 supports these functions
31 # elif defined(_MSC_VER) && _MSC_VER >= 1800
32 # define GEOGRAPHICLIB_CXX11_MATH 1
33 # else
34 # define GEOGRAPHICLIB_CXX11_MATH 0
35 # endif
36 #endif
37 
38 #if !defined(GEOGRAPHICLIB_WORDS_BIGENDIAN)
39 # define GEOGRAPHICLIB_WORDS_BIGENDIAN 0
40 #endif
41 
42 #if !defined(GEOGRAPHICLIB_HAVE_LONG_DOUBLE)
43 # define GEOGRAPHICLIB_HAVE_LONG_DOUBLE 0
44 #endif
45 
46 #if !defined(GEOGRAPHICLIB_PRECISION)
47 /**
48  * The precision of floating point numbers used in %GeographicLib. 1 means
49  * float (single precision); 2 (the default) means double; 3 means long double;
50  * 4 is reserved for quadruple precision. Nearly all the testing has been
51  * carried out with doubles and that's the recommended configuration. In order
52  * for long double to be used, GEOGRAPHICLIB_HAVE_LONG_DOUBLE needs to be
53  * defined. Note that with Microsoft Visual Studio, long double is the same as
54  * double.
55  **********************************************************************/
56 # define GEOGRAPHICLIB_PRECISION 2
57 #endif
58 
59 #include <cmath>
60 #include <algorithm>
61 #include <limits>
62 
63 #if GEOGRAPHICLIB_PRECISION == 4
64 #include <boost/multiprecision/float128.hpp>
65 #include <boost/math/special_functions/hypot.hpp>
66 #include <boost/math/special_functions/expm1.hpp>
67 #include <boost/math/special_functions/log1p.hpp>
68 #include <boost/math/special_functions/atanh.hpp>
69 #include <boost/math/special_functions/asinh.hpp>
70 #include <boost/math/special_functions/cbrt.hpp>
71 #elif GEOGRAPHICLIB_PRECISION == 5
72 #include <mpreal.h>
73 #endif
74 
75 #if GEOGRAPHICLIB_PRECISION > 3
76 // volatile keyword makes no sense for multiprec types
77 #define GEOGRAPHICLIB_VOLATILE
78 // Signal a convergence failure with multiprec types by throwing an exception
79 // at loop exit.
80 #define GEOGRAPHICLIB_PANIC \
81  (throw GeographicLib::GeographicErr("Convergence failure"), false)
82 #else
83 #define GEOGRAPHICLIB_VOLATILE volatile
84 // Ignore convergence failures with standard floating points types by allowing
85 // loop to exit cleanly.
86 #define GEOGRAPHICLIB_PANIC false
87 #endif
88 
89 namespace GeographicLib {
90 
91  /**
92  * \brief Mathematical functions needed by %GeographicLib
93  *
94  * Define mathematical functions in order to localize system dependencies and
95  * to provide generic versions of the functions. In addition define a real
96  * type to be used by %GeographicLib.
97  *
98  * Example of use:
99  * \include example-Math.cpp
100  **********************************************************************/
102  private:
103  void dummy() {
104  GEOGRAPHICLIB_STATIC_ASSERT(GEOGRAPHICLIB_PRECISION >= 1 &&
106  "Bad value of precision");
107  }
108  Math(); // Disable constructor
109  public:
110 
111 #if GEOGRAPHICLIB_HAVE_LONG_DOUBLE
112  /**
113  * The extended precision type for real numbers, used for some testing.
114  * This is long double on computers with this type; otherwise it is double.
115  **********************************************************************/
116  typedef long double extended;
117 #else
118  typedef double extended;
119 #endif
120 
121 #if GEOGRAPHICLIB_PRECISION == 2
122  /**
123  * The real type for %GeographicLib. Nearly all the testing has been done
124  * with \e real = double. However, the algorithms should also work with
125  * float and long double (where available). (<b>CAUTION</b>: reasonable
126  * accuracy typically cannot be obtained using floats.)
127  **********************************************************************/
128  typedef double real;
129 #elif GEOGRAPHICLIB_PRECISION == 1
130  typedef float real;
131 #elif GEOGRAPHICLIB_PRECISION == 3
132  typedef extended real;
133 #elif GEOGRAPHICLIB_PRECISION == 4
134  typedef boost::multiprecision::float128 real;
135 #elif GEOGRAPHICLIB_PRECISION == 5
136  typedef mpfr::mpreal real;
137 #else
138  typedef double real;
139 #endif
140 
141  /**
142  * @return the number of bits of precision in a real number.
143  **********************************************************************/
144  static inline int digits() {
145 #if GEOGRAPHICLIB_PRECISION != 5
146  return std::numeric_limits<real>::digits;
147 #else
148  return std::numeric_limits<real>::digits();
149 #endif
150  }
151 
152  /**
153  * Set the binary precision of a real number.
154  *
155  * @param[in] ndigits the number of bits of precision.
156  * @return the resulting number of bits of precision.
157  *
158  * This only has an effect when GEOGRAPHICLIB_PRECISION == 5.
159  **********************************************************************/
160  static inline int set_digits(int ndigits) {
161 #if GEOGRAPHICLIB_PRECISION != 5
162  (void)ndigits;
163 #else
164  mpfr::mpreal::set_default_prec(ndigits >= 2 ? ndigits : 2);
165 #endif
166  return digits();
167  }
168 
169  /**
170  * @return the number of decimal digits of precision in a real number.
171  **********************************************************************/
172  static inline int digits10() {
173 #if GEOGRAPHICLIB_PRECISION != 5
174  return std::numeric_limits<real>::digits10;
175 #else
176  return std::numeric_limits<real>::digits10();
177 #endif
178  }
179 
180  /**
181  * Number of additional decimal digits of precision for real relative to
182  * double (0 for float).
183  **********************************************************************/
184  static inline int extra_digits() {
185  return
186  digits10() > std::numeric_limits<double>::digits10 ?
187  digits10() - std::numeric_limits<double>::digits10 : 0;
188  }
189 
190 #if GEOGRAPHICLIB_PRECISION <= 3
191  /**
192  * Number of additional decimal digits of precision of real relative to
193  * double (0 for float).
194  *
195  * <b>DEPRECATED</b>: use extra_digits() instead
196  **********************************************************************/
197  static const int extradigits =
198  std::numeric_limits<real>::digits10 >
199  std::numeric_limits<double>::digits10 ?
200  std::numeric_limits<real>::digits10 -
201  std::numeric_limits<double>::digits10 : 0;
202 #endif
203 
204  /**
205  * true if the machine is big-endian.
206  **********************************************************************/
207  static const bool bigendian = GEOGRAPHICLIB_WORDS_BIGENDIAN;
208 
209  /**
210  * @tparam T the type of the returned value.
211  * @return &pi;.
212  **********************************************************************/
213  template<typename T> static inline T pi() {
214  using std::atan2;
215  static const T pi = atan2(T(0), T(-1));
216  return pi;
217  }
218  /**
219  * A synonym for pi<real>().
220  **********************************************************************/
221  static inline real pi() { return pi<real>(); }
222 
223  /**
224  * @tparam T the type of the returned value.
225  * @return the number of radians in a degree.
226  **********************************************************************/
227  template<typename T> static inline T degree() {
228  static const T degree = pi<T>() / 180;
229  return degree;
230  }
231  /**
232  * A synonym for degree<real>().
233  **********************************************************************/
234  static inline real degree() { return degree<real>(); }
235 
236  /**
237  * Square a number.
238  *
239  * @tparam T the type of the argument and the returned value.
240  * @param[in] x
241  * @return <i>x</i><sup>2</sup>.
242  **********************************************************************/
243  template<typename T> static inline T sq(T x)
244  { return x * x; }
245 
246  /**
247  * The hypotenuse function avoiding underflow and overflow.
248  *
249  * @tparam T the type of the arguments and the returned value.
250  * @param[in] x
251  * @param[in] y
252  * @return sqrt(<i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>).
253  **********************************************************************/
254  template<typename T> static inline T hypot(T x, T y) {
255 #if GEOGRAPHICLIB_CXX11_MATH
256  using std::hypot; return hypot(x, y);
257 #else
258  using std::abs; using std::sqrt;
259  x = abs(x); y = abs(y);
260  if (x < y) std::swap(x, y); // Now x >= y >= 0
261  y /= (x ? x : 1);
262  return x * sqrt(1 + y * y);
263  // For an alternative (square-root free) method see
264  // C. Moler and D. Morrision (1983) http://dx.doi.org/10.1147/rd.276.0577
265  // and A. A. Dubrulle (1983) http://dx.doi.org/10.1147/rd.276.0582
266 #endif
267  }
268 
269  /**
270  * exp(\e x) &minus; 1 accurate near \e x = 0.
271  *
272  * @tparam T the type of the argument and the returned value.
273  * @param[in] x
274  * @return exp(\e x) &minus; 1.
275  **********************************************************************/
276  template<typename T> static inline T expm1(T x) {
277 #if GEOGRAPHICLIB_CXX11_MATH
278  using std::expm1; return expm1(x);
279 #else
280  using std::exp; using std::abs; using std::log;
281  volatile T
282  y = exp(x),
283  z = y - 1;
284  // The reasoning here is similar to that for log1p. The expression
285  // mathematically reduces to exp(x) - 1, and the factor z/log(y) = (y -
286  // 1)/log(y) is a slowly varying quantity near y = 1 and is accurately
287  // computed.
288  return abs(x) > 1 ? z : (z == 0 ? x : x * z / log(y));
289 #endif
290  }
291 
292  /**
293  * log(1 + \e x) accurate near \e x = 0.
294  *
295  * @tparam T the type of the argument and the returned value.
296  * @param[in] x
297  * @return log(1 + \e x).
298  **********************************************************************/
299  template<typename T> static inline T log1p(T x) {
300 #if GEOGRAPHICLIB_CXX11_MATH
301  using std::log1p; return log1p(x);
302 #else
303  using std::log;
304  volatile T
305  y = 1 + x,
306  z = y - 1;
307  // Here's the explanation for this magic: y = 1 + z, exactly, and z
308  // approx x, thus log(y)/z (which is nearly constant near z = 0) returns
309  // a good approximation to the true log(1 + x)/x. The multiplication x *
310  // (log(y)/z) introduces little additional error.
311  return z == 0 ? x : x * log(y) / z;
312 #endif
313  }
314 
315  /**
316  * The inverse hyperbolic sine function.
317  *
318  * @tparam T the type of the argument and the returned value.
319  * @param[in] x
320  * @return asinh(\e x).
321  **********************************************************************/
322  template<typename T> static inline T asinh(T x) {
323 #if GEOGRAPHICLIB_CXX11_MATH
324  using std::asinh; return asinh(x);
325 #else
326  using std::abs; T y = abs(x); // Enforce odd parity
327  y = log1p(y * (1 + y/(hypot(T(1), y) + 1)));
328  return x < 0 ? -y : y;
329 #endif
330  }
331 
332  /**
333  * The inverse hyperbolic tangent function.
334  *
335  * @tparam T the type of the argument and the returned value.
336  * @param[in] x
337  * @return atanh(\e x).
338  **********************************************************************/
339  template<typename T> static inline T atanh(T x) {
340 #if GEOGRAPHICLIB_CXX11_MATH
341  using std::atanh; return atanh(x);
342 #else
343  using std::abs; T y = abs(x); // Enforce odd parity
344  y = log1p(2 * y/(1 - y))/2;
345  return x < 0 ? -y : y;
346 #endif
347  }
348 
349  /**
350  * The cube root function.
351  *
352  * @tparam T the type of the argument and the returned value.
353  * @param[in] x
354  * @return the real cube root of \e x.
355  **********************************************************************/
356  template<typename T> static inline T cbrt(T x) {
357 #if GEOGRAPHICLIB_CXX11_MATH
358  using std::cbrt; return cbrt(x);
359 #else
360  using std::abs; using std::pow;
361  T y = pow(abs(x), 1/T(3)); // Return the real cube root
362  return x < 0 ? -y : y;
363 #endif
364  }
365 
366  /**
367  * The error-free sum of two numbers.
368  *
369  * @tparam T the type of the argument and the returned value.
370  * @param[in] u
371  * @param[in] v
372  * @param[out] t the exact error given by (\e u + \e v) - \e s.
373  * @return \e s = round(\e u + \e v).
374  *
375  * See D. E. Knuth, TAOCP, Vol 2, 4.2.2, Theorem B. (Note that \e t can be
376  * the same as one of the first two arguments.)
377  **********************************************************************/
378  template<typename T> static inline T sum(T u, T v, T& t) {
379  GEOGRAPHICLIB_VOLATILE T s = u + v;
380  GEOGRAPHICLIB_VOLATILE T up = s - v;
381  GEOGRAPHICLIB_VOLATILE T vpp = s - up;
382  up -= u;
383  vpp -= v;
384  t = -(up + vpp);
385  // u + v = s + t
386  // = round(u + v) + t
387  return s;
388  }
389 
390  /**
391  * Normalize an angle (restricted input range).
392  *
393  * @tparam T the type of the argument and returned value.
394  * @param[in] x the angle in degrees.
395  * @return the angle reduced to the range [&minus;180&deg;, 180&deg;).
396  *
397  * \e x must lie in [&minus;540&deg;, 540&deg;).
398  **********************************************************************/
399  template<typename T> static inline T AngNormalize(T x)
400  { return x >= 180 ? x - 360 : (x < -180 ? x + 360 : x); }
401 
402  /**
403  * Normalize an arbitrary angle.
404  *
405  * @tparam T the type of the argument and returned value.
406  * @param[in] x the angle in degrees.
407  * @return the angle reduced to the range [&minus;180&deg;, 180&deg;).
408  *
409  * The range of \e x is unrestricted.
410  **********************************************************************/
411  template<typename T> static inline T AngNormalize2(T x)
412  { using std::fmod; return AngNormalize<T>(fmod(x, T(360))); }
413 
414  /**
415  * Difference of two angles reduced to [&minus;180&deg;, 180&deg;]
416  *
417  * @tparam T the type of the arguments and returned value.
418  * @param[in] x the first angle in degrees.
419  * @param[in] y the second angle in degrees.
420  * @return \e y &minus; \e x, reduced to the range [&minus;180&deg;,
421  * 180&deg;].
422  *
423  * \e x and \e y must both lie in [&minus;180&deg;, 180&deg;]. The result
424  * is equivalent to computing the difference exactly, reducing it to
425  * (&minus;180&deg;, 180&deg;] and rounding the result. Note that this
426  * prescription allows &minus;180&deg; to be returned (e.g., if \e x is
427  * tiny and negative and \e y = 180&deg;).
428  **********************************************************************/
429  template<typename T> static inline T AngDiff(T x, T y) {
430  T t, d = sum(-x, y, t);
431  if ((d - T(180)) + t > T(0)) // y - x > 180
432  d -= T(360); // exact
433  else if ((d + T(180)) + t <= T(0)) // y - x <= -180
434  d += T(360); // exact
435  return d + t;
436  }
437 
438  /**
439  * Test for finiteness.
440  *
441  * @tparam T the type of the argument.
442  * @param[in] x
443  * @return true if number is finite, false if NaN or infinite.
444  **********************************************************************/
445  template<typename T> static inline bool isfinite(T x) {
446 #if GEOGRAPHICLIB_CXX11_MATH
447  using std::isfinite; return isfinite(x);
448 #else
449  using std::abs;
450  return abs(x) <= (std::numeric_limits<T>::max)();
451 #endif
452  }
453 
454  /**
455  * The NaN (not a number)
456  *
457  * @tparam T the type of the returned value.
458  * @return NaN if available, otherwise return the max real of type T.
459  **********************************************************************/
460  template<typename T> static inline T NaN() {
461  return std::numeric_limits<T>::has_quiet_NaN ?
462  std::numeric_limits<T>::quiet_NaN() :
463  (std::numeric_limits<T>::max)();
464  }
465  /**
466  * A synonym for NaN<real>().
467  **********************************************************************/
468  static inline real NaN() { return NaN<real>(); }
469 
470  /**
471  * Test for NaN.
472  *
473  * @tparam T the type of the argument.
474  * @param[in] x
475  * @return true if argument is a NaN.
476  **********************************************************************/
477  template<typename T> static inline bool isnan(T x) {
478 #if GEOGRAPHICLIB_CXX11_MATH
479  using std::isnan; return isnan(x);
480 #else
481  return x != x;
482 #endif
483  }
484 
485  /**
486  * Infinity
487  *
488  * @tparam T the type of the returned value.
489  * @return infinity if available, otherwise return the max real.
490  **********************************************************************/
491  template<typename T> static inline T infinity() {
492  return std::numeric_limits<T>::has_infinity ?
493  std::numeric_limits<T>::infinity() :
494  (std::numeric_limits<T>::max)();
495  }
496  /**
497  * A synonym for infinity<real>().
498  **********************************************************************/
499  static inline real infinity() { return infinity<real>(); }
500 
501  /**
502  * Swap the bytes of a quantity
503  *
504  * @tparam T the type of the argument and the returned value.
505  * @param[in] x
506  * @return x with its bytes swapped.
507  **********************************************************************/
508  template<typename T> static inline T swab(T x) {
509  union {
510  T r;
511  unsigned char c[sizeof(T)];
512  } b;
513  b.r = x;
514  for (int i = sizeof(T)/2; i--; )
515  std::swap(b.c[i], b.c[sizeof(T) - 1 - i]);
516  return b.r;
517  }
518 
519 #if GEOGRAPHICLIB_PRECISION == 4
520  typedef boost::math::policies::policy
521  < boost::math::policies::domain_error
522  <boost::math::policies::errno_on_error>,
523  boost::math::policies::pole_error
524  <boost::math::policies::errno_on_error>,
525  boost::math::policies::overflow_error
526  <boost::math::policies::errno_on_error>,
527  boost::math::policies::evaluation_error
528  <boost::math::policies::errno_on_error> >
529  boost_special_functions_policy;
530 
531  static inline real hypot(real x, real y)
532  { return boost::math::hypot(x, y, boost_special_functions_policy()); }
533 
534  static inline real expm1(real x)
535  { return boost::math::expm1(x, boost_special_functions_policy()); }
536 
537  static inline real log1p(real x)
538  { return boost::math::log1p(x, boost_special_functions_policy()); }
539 
540  static inline real asinh(real x)
541  { return boost::math::asinh(x, boost_special_functions_policy()); }
542 
543  static inline real atanh(real x)
544  { return boost::math::atanh(x, boost_special_functions_policy()); }
545 
546  static inline real cbrt(real x)
547  { return boost::math::cbrt(x, boost_special_functions_policy()); }
548 
549  static inline bool isnan(real x) { return boost::math::isnan(x); }
550 
551  static inline bool isfinite(real x) { return boost::math::isfinite(x); }
552 #endif
553  };
554 
555 } // namespace GeographicLib
556 
557 #endif // GEOGRAPHICLIB_MATH_HPP