GeographicLib  1.21
Math.hpp
Go to the documentation of this file.
00001 /**
00002  * \file Math.hpp
00003  * \brief Header for GeographicLib::Math class
00004  *
00005  * Copyright (c) Charles Karney (2008-2011) <charles@karney.com> and licensed
00006  * under the MIT/X11 License.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 // Constants.hpp includes Math.hpp.  Place this include outside Math.hpp's
00011 // include guard to enforce this ordering.
00012 #include <GeographicLib/Constants.hpp>
00013 
00014 #if !defined(GEOGRAPHICLIB_MATH_HPP)
00015 #define GEOGRAPHICLIB_MATH_HPP "$Id: edd244e4c5c74e696096c2b6d598728957a0d36d $"
00016 
00017 /**
00018  * Are C++11 math functions available?
00019  **********************************************************************/
00020 #if !defined(GEOGRAPHICLIB_CPLUSPLUS11_MATH)
00021 #  if defined(__GXX_EXPERIMENTAL_CXX0X__)
00022 #    define GEOGRAPHICLIB_CPLUSPLUS11_MATH 1
00023 #  else
00024 #    define GEOGRAPHICLIB_CPLUSPLUS11_MATH 0
00025 #  endif
00026 #endif
00027 
00028 #if !defined(WORDS_BIGENDIAN)
00029 # define WORDS_BIGENDIAN 0
00030 #endif
00031 
00032 #if !defined(GEOGRAPHICLIB_PREC)
00033 /**
00034  * The precision of floating point numbers used in %GeographicLib.  0 means
00035  * float; 1 (default) means double; 2 means long double.  Nearly all the
00036  * testing has been carried out with doubles and that's the recommended
00037  * configuration.  In order for long double to be used, HAVE_LONG_DOUBLE needs
00038  * to be defined.  Note that with Microsoft Visual Studio, long double is the
00039  * same as double.
00040  **********************************************************************/
00041 #define GEOGRAPHICLIB_PREC 1
00042 #endif
00043 
00044 #include <cmath>
00045 #include <limits>
00046 #include <algorithm>
00047 #include <vector>
00048 
00049 namespace GeographicLib {
00050 
00051   /**
00052    * \brief Mathematical functions needed by %GeographicLib
00053    *
00054    * Define mathematical functions in order to localize system dependencies and
00055    * to provide generic versions of the functions.  In addition define a real
00056    * type to be used by %GeographicLib.
00057    *
00058    * Example of use:
00059    * \include example-Math.cpp
00060    **********************************************************************/
00061   class GEOGRAPHIC_EXPORT Math {
00062   private:
00063     void dummy() {
00064       STATIC_ASSERT(GEOGRAPHICLIB_PREC >= 0 && GEOGRAPHICLIB_PREC <= 2,
00065                     "Bad value of precision");
00066     }
00067     Math();                     // Disable constructor
00068   public:
00069 
00070 #if defined(HAVE_LONG_DOUBLE)
00071     /**
00072      * The extended precision type for real numbers, used for some testing.
00073      * This is long double on computers with this type; otherwise it is double.
00074      **********************************************************************/
00075     typedef long double extended;
00076 #else
00077     typedef double extended;
00078 #endif
00079 
00080 #if GEOGRAPHICLIB_PREC == 1
00081     /**
00082      * The real type for %GeographicLib. Nearly all the testing has been done
00083      * with \e real = double.  However, the algorithms should also work with
00084      * float and long double (where available).  (<b>CAUTION</b>: reasonable
00085      * accuracy typically cannot be obtained using floats.)
00086      **********************************************************************/
00087     typedef double real;
00088 #elif GEOGRAPHICLIB_PREC == 0
00089     typedef float real;
00090 #elif GEOGRAPHICLIB_PREC == 2
00091     typedef extended real;
00092 #else
00093     typedef double real;
00094 #endif
00095 
00096     /**
00097      * true if the machine is big-endian
00098      **********************************************************************/
00099     static const bool bigendian = WORDS_BIGENDIAN;
00100 
00101     /**
00102      * @tparam T the type of the returned value.
00103      * @return \e pi.
00104      **********************************************************************/
00105     template<typename T> static inline T pi() throw()
00106     { return std::atan2(T(0), -T(1)); }
00107     /**
00108      * A synonym for pi<real>().
00109      **********************************************************************/
00110     static inline real pi() throw() { return pi<real>(); }
00111 
00112     /**
00113      * @tparam T the type of the returned value.
00114      * @return the number of radians in a degree.
00115      **********************************************************************/
00116     template<typename T> static inline T degree() throw()
00117     { return pi<T>() / T(180); }
00118     /**
00119      * A synonym for degree<real>().
00120      **********************************************************************/
00121     static inline real degree() throw() { return degree<real>(); }
00122 
00123     /**
00124      * Square a number.
00125 
00126      * @tparam T the type of the argument and the returned value.
00127      * @param[in] x
00128      * @return <i>x</i><sup>2</sup>.
00129      **********************************************************************/
00130     template<typename T> static inline T sq(T x) throw()
00131     { return x * x; }
00132 
00133 #if defined(DOXYGEN)
00134     /**
00135      * The hypotenuse function avoiding underflow and overflow.
00136      *
00137      * @tparam T the type of the arguments and the returned value.
00138      * @param[in] x
00139      * @param[in] y
00140      * @return sqrt(<i>x</i><sup>2</sup> + <i>y</i><sup>2</sup>).
00141      **********************************************************************/
00142     template<typename T> static inline T hypot(T x, T y) throw() {
00143       x = std::abs(x);
00144       y = std::abs(y);
00145       T a = (std::max)(x, y),
00146         b = (std::min)(x, y) / (a ? a : 1);
00147       return a * std::sqrt(1 + b * b);
00148     }
00149 #elif GEOGRAPHICLIB_CPLUSPLUS11_MATH
00150     template<typename T> static inline T hypot(T x, T y) throw()
00151     { return std::hypot(x, y); }
00152 #elif defined(_MSC_VER)
00153     static inline double hypot(double x, double y) throw()
00154     { return _hypot(x, y); }
00155 #if _MSC_VER < 1400
00156     // Visual C++ 7.1/VS .NET 2003 does not have _hypotf()
00157     static inline float hypot(float x, float y) throw()
00158     { return float(_hypot(x, y)); }
00159 #else
00160     static inline float hypot(float x, float y) throw()
00161     { return _hypotf(x, y); }
00162 #endif
00163 #if defined(HAVE_LONG_DOUBLE)
00164     static inline long double hypot(long double x, long double y) throw()
00165     { return _hypot(x, y); }
00166 #endif
00167 #else
00168     // Use overloading to define generic versions
00169     static inline double hypot(double x, double y) throw()
00170     { return ::hypot(x, y); }
00171     static inline float hypot(float x, float y) throw()
00172     { return ::hypotf(x, y); }
00173 #if defined(HAVE_LONG_DOUBLE)
00174     static inline long double hypot(long double x, long double y) throw()
00175     { return ::hypotl(x, y); }
00176 #endif
00177 #endif
00178 
00179 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
00180     /**
00181      * exp(\e x) - 1 accurate near \e x = 0.  This is taken from
00182      * N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd
00183      * Edition (SIAM, 2002), Sec 1.14.1, p 19.
00184      *
00185      * @tparam T the type of the argument and the returned value.
00186      * @param[in] x
00187      * @return exp(\e x) - 1.
00188      **********************************************************************/
00189     template<typename T> static inline T expm1(T x) throw() {
00190       volatile T
00191         y = std::exp(x),
00192         z = y - 1;
00193       // The reasoning here is similar to that for log1p.  The expression
00194       // mathematically reduces to exp(x) - 1, and the factor z/log(y) = (y -
00195       // 1)/log(y) is a slowly varying quantity near y = 1 and is accurately
00196       // computed.
00197       return std::abs(x) > 1 ? z : (z == 0 ? x : x * z / std::log(y));
00198     }
00199 #elif GEOGRAPHICLIB_CPLUSPLUS11_MATH
00200     template<typename T> static inline T expm1(T x) throw()
00201     { return std::expm1(x); }
00202 #else
00203     static inline double expm1(double x) throw() { return ::expm1(x); }
00204     static inline float expm1(float x) throw() { return ::expm1f(x); }
00205 #if defined(HAVE_LONG_DOUBLE)
00206     static inline long double expm1(long double x) throw()
00207     { return ::expm1l(x); }
00208 #endif
00209 #endif
00210 
00211 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
00212     /**
00213      * log(1 + \e x) accurate near \e x = 0.
00214      *
00215      * This is taken from D. Goldberg,
00216      * <a href="http://dx.doi.org/10.1145/103162.103163">What every computer
00217      * scientist should know about floating-point arithmetic</a> (1991),
00218      * Theorem 4.  See also, Higham (op. cit.), Answer to Problem 1.5, p 528.
00219      *
00220      * @tparam T the type of the argument and the returned value.
00221      * @param[in] x
00222      * @return log(1 + \e x).
00223      **********************************************************************/
00224     template<typename T> static inline T log1p(T x) throw() {
00225       volatile T
00226         y = 1 + x,
00227         z = y - 1;
00228       // Here's the explanation for this magic: y = 1 + z, exactly, and z
00229       // approx x, thus log(y)/z (which is nearly constant near z = 0) returns
00230       // a good approximation to the true log(1 + x)/x.  The multiplication x *
00231       // (log(y)/z) introduces little additional error.
00232       return z == 0 ? x : x * std::log(y) / z;
00233     }
00234 #elif GEOGRAPHICLIB_CPLUSPLUS11_MATH
00235     template<typename T> static inline T log1p(T x) throw()
00236     { return std::log1p(x); }
00237 #else
00238     static inline double log1p(double x) throw() { return ::log1p(x); }
00239     static inline float log1p(float x) throw() { return ::log1pf(x); }
00240 #if defined(HAVE_LONG_DOUBLE)
00241     static inline long double log1p(long double x) throw()
00242     { return ::log1pl(x); }
00243 #endif
00244 #endif
00245 
00246 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
00247     /**
00248      * The inverse hyperbolic sine function.  This is defined in terms of
00249      * Math::log1p(\e x) in order to maintain accuracy near \e x = 0.  In
00250      * addition, the odd parity of the function is enforced.
00251      *
00252      * @tparam T the type of the argument and the returned value.
00253      * @param[in] x
00254      * @return asinh(\e x).
00255      **********************************************************************/
00256     template<typename T> static inline T asinh(T x) throw() {
00257       T y = std::abs(x);     // Enforce odd parity
00258       y = log1p(y * (1 + y/(hypot(T(1), y) + 1)));
00259       return x < 0 ? -y : y;
00260     }
00261 #elif GEOGRAPHICLIB_CPLUSPLUS11_MATH
00262     template<typename T> static inline T asinh(T x) throw()
00263     { return std::asinh(x); }
00264 #else
00265     static inline double asinh(double x) throw() { return ::asinh(x); }
00266     static inline float asinh(float x) throw() { return ::asinhf(x); }
00267 #if defined(HAVE_LONG_DOUBLE)
00268     static inline long double asinh(long double x) throw()
00269     { return ::asinhl(x); }
00270 #endif
00271 #endif
00272 
00273 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
00274     /**
00275      * The inverse hyperbolic tangent function.  This is defined in terms of
00276      * Math::log1p(\e x) in order to maintain accuracy near \e x = 0.  In
00277      * addition, the odd parity of the function is enforced.
00278      *
00279      * @tparam T the type of the argument and the returned value.
00280      * @param[in] x
00281      * @return atanh(\e x).
00282      **********************************************************************/
00283     template<typename T> static inline T atanh(T x) throw() {
00284       T y = std::abs(x);     // Enforce odd parity
00285       y = log1p(2 * y/(1 - y))/2;
00286       return x < 0 ? -y : y;
00287     }
00288 #elif GEOGRAPHICLIB_CPLUSPLUS11_MATH
00289     template<typename T> static inline T atanh(T x) throw()
00290     { return std::atanh(x); }
00291 #else
00292     static inline double atanh(double x) throw() { return ::atanh(x); }
00293     static inline float atanh(float x) throw() { return ::atanhf(x); }
00294 #if defined(HAVE_LONG_DOUBLE)
00295     static inline long double atanh(long double x) throw()
00296     { return ::atanhl(x); }
00297 #endif
00298 #endif
00299 
00300 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
00301     /**
00302      * The cube root function.
00303      *
00304      * @tparam T the type of the argument and the returned value.
00305      * @param[in] x
00306      * @return the real cube root of \e x.
00307      **********************************************************************/
00308     template<typename T> static inline T cbrt(T x) throw() {
00309       T y = std::pow(std::abs(x), 1/T(3)); // Return the real cube root
00310       return x < 0 ? -y : y;
00311     }
00312 #elif GEOGRAPHICLIB_CPLUSPLUS11_MATH
00313     template<typename T> static inline T cbrt(T x) throw()
00314     { return std::cbrt(x); }
00315 #else
00316     static inline double cbrt(double x) throw() { return ::cbrt(x); }
00317     static inline float cbrt(float x) throw() { return ::cbrtf(x); }
00318 #if defined(HAVE_LONG_DOUBLE)
00319     static inline long double cbrt(long double x) throw() { return ::cbrtl(x); }
00320 #endif
00321 #endif
00322 
00323     /**
00324      * Test for finiteness.
00325      *
00326      * @tparam T the type of the argument.
00327      * @param[in] x
00328      * @return true if number is finite, false if NaN or infinite.
00329      **********************************************************************/
00330     template<typename T> static inline bool isfinite(T x) throw() {
00331 #if defined(DOXYGEN)
00332       return std::abs(x) <= (std::numeric_limits<T>::max)();
00333 #elif (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
00334       return _finite(x) != 0;
00335 #else
00336       return std::isfinite(x);
00337 #endif
00338     }
00339 
00340     /**
00341      * The NaN (not a number)
00342      *
00343      * @tparam T the type of the returned value.
00344      * @return NaN if available, otherwise return the max real.
00345      **********************************************************************/
00346     template<typename T> static inline T NaN() throw() {
00347       return std::numeric_limits<T>::has_quiet_NaN ?
00348         std::numeric_limits<T>::quiet_NaN() :
00349         (std::numeric_limits<T>::max)();
00350     }
00351     /**
00352      * A synonym for NaN<real>().
00353      **********************************************************************/
00354     static inline real NaN() throw() { return NaN<real>(); }
00355 
00356     /**
00357      * Test for NaN.
00358      *
00359      * @tparam T the type of the argument.
00360      * @param[in] x
00361      * @return true if argument is a NaN.
00362      **********************************************************************/
00363     template<typename T> static inline bool isnan(T x) throw() {
00364 #if defined(DOXYGEN) || (defined(_MSC_VER) && !GEOGRAPHICLIB_CPLUSPLUS11_MATH)
00365       return x != x;
00366 #else
00367       return std::isnan(x);
00368 #endif
00369     }
00370 
00371     /**
00372      * Infinity
00373      *
00374      * @tparam T the type of the returned value.
00375      * @return infinity if available, otherwise return the max real.
00376      **********************************************************************/
00377     template<typename T> static inline T infinity() throw() {
00378       return std::numeric_limits<T>::has_infinity ?
00379         std::numeric_limits<T>::infinity() :
00380         (std::numeric_limits<T>::max)();
00381     }
00382     /**
00383      * A synonym for infinity<real>().
00384      **********************************************************************/
00385     static inline real infinity() throw() { return infinity<real>(); }
00386 
00387     /**
00388      * Swap the bytes of a quantity
00389      *
00390      * @tparam T the type of the argument and the returned value.
00391      * @param[in] x
00392      * @return x with its bytes swapped.
00393      **********************************************************************/
00394     template<typename T> static inline T swab(T x) {
00395       union {
00396         T r;
00397         unsigned char c[sizeof(T)];
00398       } b;
00399       b.r = x;
00400       for (int i = sizeof(T)/2; i--; )
00401         std::swap(b.c[i], b.c[sizeof(T) - 1 - i]);
00402       return b.r;
00403     }
00404   };
00405 
00406 } // namespace GeographicLib
00407 
00408 #endif  // GEOGRAPHICLIB_MATH_HPP