GeographicLib  1.35
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-2011) <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_CPLUSPLUS11_MATH)
21 # if defined(__GNUC__) && __GNUC__ == 4 && __GNUC_MINOR__ >= 7 \
22  && __cplusplus >= 201103 && !(defined(__ANDROID__) || defined(ANDROID))
23 // The android toolchain uses g++ and supports C++11, but not, apparently, the
24 // new mathematical functions introduced with C++11. Android toolchains might
25 // define __ANDROID__ or ANDROID; so need to check both.
26 # define GEOGRAPHICLIB_CPLUSPLUS11_MATH 1
27 # elif defined(_MSC_VER) && _MSC_VER >= 1800
28 # define GEOGRAPHICLIB_CPLUSPLUS11_MATH 1
29 # else
30 # define GEOGRAPHICLIB_CPLUSPLUS11_MATH 0
31 # endif
32 #endif
33 
34 #if !defined(WORDS_BIGENDIAN)
35 # define WORDS_BIGENDIAN 0
36 #endif
37 
38 #if !defined(HAVE_LONG_DOUBLE)
39 # define HAVE_LONG_DOUBLE 0
40 #endif
41 
42 #if !defined(GEOGRAPHICLIB_PRECISION)
43 /**
44  * The precision of floating point numbers used in %GeographicLib. 1 means
45  * float (single precision); 2 (the default) means double; 3 means long double;
46  * 4 is reserved for quadruple precision. Nearly all the testing has been
47  * carried out with doubles and that's the recommended configuration. In order
48  * for long double to be used, HAVE_LONG_DOUBLE needs to be defined. Note that
49  * with Microsoft Visual Studio, long double is the same as double.
50  **********************************************************************/
51 # define GEOGRAPHICLIB_PRECISION 2
52 #endif
53 
54 #include <cmath>
55 #include <algorithm>
56 #include <limits>
57 #if defined(_LIBCPP_VERSION)
58 #include <type_traits>
59 #endif
60 
61 namespace GeographicLib {
62 
63  /**
64  * \brief Mathematical functions needed by %GeographicLib
65  *
66  * Define mathematical functions in order to localize system dependencies and
67  * to provide generic versions of the functions. In addition define a real
68  * type to be used by %GeographicLib.
69  *
70  * Example of use:
71  * \include example-Math.cpp
72  **********************************************************************/
74  private:
75  void dummy() {
78  "Bad value of precision");
79  }
80  Math(); // Disable constructor
81  public:
82 
83 #if HAVE_LONG_DOUBLE
84  /**
85  * The extended precision type for real numbers, used for some testing.
86  * This is long double on computers with this type; otherwise it is double.
87  **********************************************************************/
88  typedef long double extended;
89 #else
90  typedef double extended;
91 #endif
92 
93 #if GEOGRAPHICLIB_PRECISION == 2
94  /**
95  * The real type for %GeographicLib. Nearly all the testing has been done
96  * with \e real = double. However, the algorithms should also work with
97  * float and long double (where available). (<b>CAUTION</b>: reasonable
98  * accuracy typically cannot be obtained using floats.)
99  **********************************************************************/
100  typedef double real;
101 #elif GEOGRAPHICLIB_PRECISION == 1
102  typedef float real;
103 #elif GEOGRAPHICLIB_PRECISION == 3
104  typedef extended real;
105 #else
106  typedef double real;
107 #endif
108 
109  /**
110  * Number of additional decimal digits of precision of real relative to
111  * double (0 for float).
112  **********************************************************************/
113  static const int extradigits =
114  std::numeric_limits<real>::digits10 >
115  std::numeric_limits<double>::digits10 ?
116  std::numeric_limits<real>::digits10 -
117  std::numeric_limits<double>::digits10 : 0;
118 
119  /**
120  * true if the machine is big-endian.
121  **********************************************************************/
122  static const bool bigendian = WORDS_BIGENDIAN;
123 
124  /**
125  * @tparam T the type of the returned value.
126  * @return &pi;.
127  **********************************************************************/
128  template<typename T> static inline T pi() throw()
129  { return std::atan2(T(0), -T(1)); }
130  /**
131  * A synonym for pi<real>().
132  **********************************************************************/
133  static inline real pi() throw() { return pi<real>(); }
134 
135  /**
136  * @tparam T the type of the returned value.
137  * @return the number of radians in a degree.
138  **********************************************************************/
139  template<typename T> static inline T degree() throw()
140  { return pi<T>() / T(180); }
141  /**
142  * A synonym for degree<real>().
143  **********************************************************************/
144  static inline real degree() throw() { return degree<real>(); }
145 
146  /**
147  * Square a number.
148  *
149  * @tparam T the type of the argument and the returned value.
150  * @param[in] x
151  * @return <i>x</i><sup>2</sup>.
152  **********************************************************************/
153  template<typename T> static inline T sq(T x) throw()
154  { return x * x; }
155 
156 #if defined(DOXYGEN)
157  /**
158  * The hypotenuse function avoiding underflow and overflow.
159  *
160  * @tparam T the type of the arguments and the returned value.
161  * @param[in] x
162  * @param[in] y
163  * @return sqrt(<i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>).
164  **********************************************************************/
165  template<typename T> static inline T hypot(T x, T y) throw() {
166  x = std::abs(x); y = std::abs(y);
167  T a = (std::max)(x, y), b = (std::min)(x, y) / (a ? a : 1);
168  return a * std::sqrt(1 + b * b);
169  // For an alternative (square-root free) method see
170  // C. Moler and D. Morrision (1983) http://dx.doi.org/10.1147/rd.276.0577
171  // and A. A. Dubrulle (1983) http://dx.doi.org/10.1147/rd.276.0582
172  }
173 #elif GEOGRAPHICLIB_CPLUSPLUS11_MATH || (defined(_MSC_VER) && _MSC_VER >= 1700)
174  template<typename T> static inline T hypot(T x, T y) throw()
175  { return std::hypot(x, y); }
176 # if HAVE_LONG_DOUBLE && defined(_MSC_VER) && _MSC_VER == 1700
177  // Visual C++ 11 doesn't have a long double overload for std::hypot --
178  // reported to MS on 2013-07-18
179  // http://connect.microsoft.com/VisualStudio/feedback/details/794416
180  // suppress the resulting "loss of data warning" with
181  static inline long double hypot(long double x, long double y) throw()
182  { return std::hypot(double(x), double(y)); }
183 # endif
184 #elif defined(_MSC_VER)
185  static inline double hypot(double x, double y) throw()
186  { return _hypot(x, y); }
187 # if _MSC_VER < 1400
188  // Visual C++ 7.1/VS .NET 2003 does not have _hypotf()
189  static inline float hypot(float x, float y) throw()
190  { return float(_hypot(x, y)); }
191 # else
192  static inline float hypot(float x, float y) throw()
193  { return _hypotf(x, y); }
194 # endif
195 # if HAVE_LONG_DOUBLE
196  static inline long double hypot(long double x, long double y) throw()
197  { return _hypot(double(x), double(y)); } // Suppress loss of data warning
198 # endif
199 #else
200  // Use overloading to define generic versions
201  static inline double hypot(double x, double y) throw()
202  { return ::hypot(x, y); }
203  static inline float hypot(float x, float y) throw()
204  { return ::hypotf(x, y); }
205 # if HAVE_LONG_DOUBLE
206  static inline long double hypot(long double x, long double y) throw()
207  { return ::hypotl(x, y); }
208 # endif
209 #endif
210 
211 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
212  /**
213  * exp(\e x) &minus; 1 accurate near \e x = 0. This is taken from
214  * N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd
215  * Edition (SIAM, 2002), Sec 1.14.1, p 19.
216  *
217  * @tparam T the type of the argument and the returned value.
218  * @param[in] x
219  * @return exp(\e x) &minus; 1.
220  **********************************************************************/
221  template<typename T> static inline T expm1(T x) throw() {
222  volatile T
223  y = std::exp(x),
224  z = y - 1;
225  // The reasoning here is similar to that for log1p. The expression
226  // mathematically reduces to exp(x) - 1, and the factor z/log(y) = (y -
227  // 1)/log(y) is a slowly varying quantity near y = 1 and is accurately
228  // computed.
229  return std::abs(x) > 1 ? z : (z == 0 ? x : x * z / std::log(y));
230  }
231 #elif GEOGRAPHICLIB_CPLUSPLUS11_MATH
232  template<typename T> static inline T expm1(T x) throw()
233  { return std::expm1(x); }
234 #else
235  static inline double expm1(double x) throw() { return ::expm1(x); }
236  static inline float expm1(float x) throw() { return ::expm1f(x); }
237 # if HAVE_LONG_DOUBLE
238  static inline long double expm1(long double x) throw()
239  { return ::expm1l(x); }
240 # endif
241 #endif
242 
243 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
244  /**
245  * log(1 + \e x) accurate near \e x = 0.
246  *
247  * This is taken from D. Goldberg,
248  * <a href="http://dx.doi.org/10.1145/103162.103163">What every computer
249  * scientist should know about floating-point arithmetic</a> (1991),
250  * Theorem 4. See also, Higham (op. cit.), Answer to Problem 1.5, p 528.
251  *
252  * @tparam T the type of the argument and the returned value.
253  * @param[in] x
254  * @return log(1 + \e x).
255  **********************************************************************/
256  template<typename T> static inline T log1p(T x) throw() {
257  volatile T
258  y = 1 + x,
259  z = y - 1;
260  // Here's the explanation for this magic: y = 1 + z, exactly, and z
261  // approx x, thus log(y)/z (which is nearly constant near z = 0) returns
262  // a good approximation to the true log(1 + x)/x. The multiplication x *
263  // (log(y)/z) introduces little additional error.
264  return z == 0 ? x : x * std::log(y) / z;
265  }
266 #elif GEOGRAPHICLIB_CPLUSPLUS11_MATH
267  template<typename T> static inline T log1p(T x) throw()
268  { return std::log1p(x); }
269 #else
270  static inline double log1p(double x) throw() { return ::log1p(x); }
271  static inline float log1p(float x) throw() { return ::log1pf(x); }
272 # if HAVE_LONG_DOUBLE
273  static inline long double log1p(long double x) throw()
274  { return ::log1pl(x); }
275 # endif
276 #endif
277 
278 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
279  /**
280  * The inverse hyperbolic sine function. This is defined in terms of
281  * Math::log1p(\e x) in order to maintain accuracy near \e x = 0. In
282  * addition, the odd parity of the function is enforced.
283  *
284  * @tparam T the type of the argument and the returned value.
285  * @param[in] x
286  * @return asinh(\e x).
287  **********************************************************************/
288  template<typename T> static inline T asinh(T x) throw() {
289  T y = std::abs(x); // Enforce odd parity
290  y = log1p(y * (1 + y/(hypot(T(1), y) + 1)));
291  return x < 0 ? -y : y;
292  }
293 #elif GEOGRAPHICLIB_CPLUSPLUS11_MATH
294  template<typename T> static inline T asinh(T x) throw()
295  { return std::asinh(x); }
296 #else
297  static inline double asinh(double x) throw() { return ::asinh(x); }
298  static inline float asinh(float x) throw() { return ::asinhf(x); }
299 # if HAVE_LONG_DOUBLE
300  static inline long double asinh(long double x) throw()
301  { return ::asinhl(x); }
302 # endif
303 #endif
304 
305 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
306  /**
307  * The inverse hyperbolic tangent function. This is defined in terms of
308  * Math::log1p(\e x) in order to maintain accuracy near \e x = 0. In
309  * addition, the odd parity of the function is enforced.
310  *
311  * @tparam T the type of the argument and the returned value.
312  * @param[in] x
313  * @return atanh(\e x).
314  **********************************************************************/
315  template<typename T> static inline T atanh(T x) throw() {
316  T y = std::abs(x); // Enforce odd parity
317  y = log1p(2 * y/(1 - y))/2;
318  return x < 0 ? -y : y;
319  }
320 #elif GEOGRAPHICLIB_CPLUSPLUS11_MATH
321  template<typename T> static inline T atanh(T x) throw()
322  { return std::atanh(x); }
323 #else
324  static inline double atanh(double x) throw() { return ::atanh(x); }
325  static inline float atanh(float x) throw() { return ::atanhf(x); }
326 # if HAVE_LONG_DOUBLE
327  static inline long double atanh(long double x) throw()
328  { return ::atanhl(x); }
329 # endif
330 #endif
331 
332 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
333  /**
334  * The cube root function.
335  *
336  * @tparam T the type of the argument and the returned value.
337  * @param[in] x
338  * @return the real cube root of \e x.
339  **********************************************************************/
340  template<typename T> static inline T cbrt(T x) throw() {
341  T y = std::pow(std::abs(x), 1/T(3)); // Return the real cube root
342  return x < 0 ? -y : y;
343  }
344 #elif GEOGRAPHICLIB_CPLUSPLUS11_MATH
345  template<typename T> static inline T cbrt(T x) throw()
346  { return std::cbrt(x); }
347 #else
348  static inline double cbrt(double x) throw() { return ::cbrt(x); }
349  static inline float cbrt(float x) throw() { return ::cbrtf(x); }
350 # if HAVE_LONG_DOUBLE
351  static inline long double cbrt(long double x) throw() { return ::cbrtl(x); }
352 # endif
353 #endif
354 
355  /**
356  * The error-free sum of two numbers.
357  *
358  * @tparam T the type of the argument and the returned value.
359  * @param[in] u
360  * @param[in] v
361  * @param[out] t the exact error given by (\e u + \e v) - \e s.
362  * @return \e s = round(\e u + \e v).
363  *
364  * See D. E. Knuth, TAOCP, Vol 2, 4.2.2, Theorem B. (Note that \e t can be
365  * the same as one of the first two arguments.)
366  **********************************************************************/
367  template<typename T> static inline T sum(T u, T v, T& t) throw() {
368  volatile T s = u + v;
369  volatile T up = s - v;
370  volatile T vpp = s - up;
371  up -= u;
372  vpp -= v;
373  t = -(up + vpp);
374  // u + v = s + t
375  // = round(u + v) + t
376  return s;
377  }
378 
379  /**
380  * Normalize an angle (restricted input range).
381  *
382  * @tparam T the type of the argument and returned value.
383  * @param[in] x the angle in degrees.
384  * @return the angle reduced to the range [&minus;180&deg;, 180&deg;).
385  *
386  * \e x must lie in [&minus;540&deg;, 540&deg;).
387  **********************************************************************/
388  template<typename T> static inline T AngNormalize(T x) throw()
389  { return x >= 180 ? x - 360 : (x < -180 ? x + 360 : x); }
390 
391  /**
392  * Normalize an arbitrary angle.
393  *
394  * @tparam T the type of the argument and returned value.
395  * @param[in] x the angle in degrees.
396  * @return the angle reduced to the range [&minus;180&deg;, 180&deg;).
397  *
398  * The range of \e x is unrestricted.
399  **********************************************************************/
400  template<typename T> static inline T AngNormalize2(T x) throw()
401  { return AngNormalize<T>(std::fmod(x, T(360))); }
402 
403  /**
404  * Difference of two angles reduced to [&minus;180&deg;, 180&deg;]
405  *
406  * @tparam T the type of the arguments and returned value.
407  * @param[in] x the first angle in degrees.
408  * @param[in] y the second angle in degrees.
409  * @return \e y &minus; \e x, reduced to the range [&minus;180&deg;,
410  * 180&deg;].
411  *
412  * \e x and \e y must both lie in [&minus;180&deg;, 180&deg;]. The result
413  * is equivalent to computing the difference exactly, reducing it to
414  * (&minus;180&deg;, 180&deg;] and rounding the result. Note that this
415  * prescription allows &minus;180&deg; to be returned (e.g., if \e x is
416  * tiny and negative and \e y = 180&deg;).
417  **********************************************************************/
418  template<typename T> static inline T AngDiff(T x, T y) throw() {
419  T t, d = sum(-x, y, t);
420  if ((d - T(180)) + t > T(0)) // y - x > 180
421  d -= T(360); // exact
422  else if ((d + T(180)) + t <= T(0)) // y - x <= -180
423  d += T(360); // exact
424  return d + t;
425  }
426 
427 #if defined(DOXYGEN)
428  /**
429  * Test for finiteness.
430  *
431  * @tparam T the type of the argument.
432  * @param[in] x
433  * @return true if number is finite, false if NaN or infinite.
434  **********************************************************************/
435  template<typename T> static inline bool isfinite(T x) throw() {
436  return std::abs(x) <= (std::numeric_limits<T>::max)();
437  }
438 #elif (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
439  template<typename T> static inline bool isfinite(T x) throw() {
440  return _finite(double(x)) != 0;
441  }
442 #elif defined(_LIBCPP_VERSION)
443  // libc++ implements std::isfinite() as a template that only allows
444  // floating-point types. isfinite is invoked by Utility::str to format
445  // numbers conveniently and this allows integer arguments, so we need to
446  // allow Math::isfinite to work on integers.
447  template<typename T> static inline
448  typename std::enable_if<std::is_floating_point<T>::value, bool>::type
449  isfinite(T x) throw() {
450  return std::isfinite(x);
451  }
452  template<typename T> static inline
453  typename std::enable_if<!std::is_floating_point<T>::value, bool>::type
454  isfinite(T /*x*/) throw() {
455  return true;
456  }
457 #else
458  template<typename T> static inline bool isfinite(T x) throw() {
459  return std::isfinite(x);
460  }
461 #endif
462 
463  /**
464  * The NaN (not a number)
465  *
466  * @tparam T the type of the returned value.
467  * @return NaN if available, otherwise return the max real of type T.
468  **********************************************************************/
469  template<typename T> static inline T NaN() throw() {
470  return std::numeric_limits<T>::has_quiet_NaN ?
471  std::numeric_limits<T>::quiet_NaN() :
472  (std::numeric_limits<T>::max)();
473  }
474  /**
475  * A synonym for NaN<real>().
476  **********************************************************************/
477  static inline real NaN() throw() { return NaN<real>(); }
478 
479  /**
480  * Test for NaN.
481  *
482  * @tparam T the type of the argument.
483  * @param[in] x
484  * @return true if argument is a NaN.
485  **********************************************************************/
486  template<typename T> static inline bool isnan(T x) throw() {
487 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
488  return x != x;
489 #else
490  return std::isnan(x);
491 #endif
492  }
493 
494  /**
495  * Infinity
496  *
497  * @tparam T the type of the returned value.
498  * @return infinity if available, otherwise return the max real.
499  **********************************************************************/
500  template<typename T> static inline T infinity() throw() {
501  return std::numeric_limits<T>::has_infinity ?
502  std::numeric_limits<T>::infinity() :
503  (std::numeric_limits<T>::max)();
504  }
505  /**
506  * A synonym for infinity<real>().
507  **********************************************************************/
508  static inline real infinity() throw() { return infinity<real>(); }
509 
510  /**
511  * Swap the bytes of a quantity
512  *
513  * @tparam T the type of the argument and the returned value.
514  * @param[in] x
515  * @return x with its bytes swapped.
516  **********************************************************************/
517  template<typename T> static inline T swab(T x) {
518  union {
519  T r;
520  unsigned char c[sizeof(T)];
521  } b;
522  b.r = x;
523  for (int i = sizeof(T)/2; i--; )
524  std::swap(b.c[i], b.c[sizeof(T) - 1 - i]);
525  return b.r;
526  }
527  };
528 
529 } // namespace GeographicLib
530 
531 #endif // GEOGRAPHICLIB_MATH_HPP
static T AngNormalize(T x)
Definition: Math.hpp:388
static T NaN()
Definition: Math.hpp:469
static T sum(T u, T v, T &t)
Definition: Math.hpp:367
static T pi()
Definition: Math.hpp:128
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:52
static T infinity()
Definition: Math.hpp:500
static real degree()
Definition: Math.hpp:144
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
static T cbrt(T x)
Definition: Math.hpp:340
static bool isfinite(T x)
Definition: Math.hpp:435
static bool isnan(T x)
Definition: Math.hpp:486
#define WORDS_BIGENDIAN
Definition: Math.hpp:35
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:73
double extended
Definition: Math.hpp:90
static T atanh(T x)
Definition: Math.hpp:315
static T expm1(T x)
Definition: Math.hpp:221
#define GEOGRAPHICLIB_PRECISION
Definition: Math.hpp:51
static T asinh(T x)
Definition: Math.hpp:288
static T hypot(T x, T y)
Definition: Math.hpp:165
static T sq(T x)
Definition: Math.hpp:153
static real pi()
Definition: Math.hpp:133
static T degree()
Definition: Math.hpp:139
static T AngDiff(T x, T y)
Definition: Math.hpp:418
#define STATIC_ASSERT(cond, reason)
Definition: Constants.hpp:37
static real NaN()
Definition: Math.hpp:477
static T log1p(T x)
Definition: Math.hpp:256
static T swab(T x)
Definition: Math.hpp:517
static real infinity()
Definition: Math.hpp:508
Header for GeographicLib::Constants class.
static T AngNormalize2(T x)
Definition: Math.hpp:400