root/Math/Math.h

0001 // @(#)root/mathcore:$Id$
0002 // Author: L. Moneta Tue Nov 14 15:44:38 2006
0003
0004 /**********************************************************************
0005  *                                                                    *
0006  * Copyright (c) 2006  LCG ROOT Math Team, CERN/PH-SFT                *
0007  *                                                                    *
0008  *                                                                    *
0009  **********************************************************************/
0010
0011 // mathematical constants like Pi
0012
0013 #ifndef ROOT_Math_Math
0014 #define ROOT_Math_Math
0015
0016 #ifdef _MSC_VER
0017 #define _USE_MATH_DEFINES
0018 #endif
0019
0020 #include <cmath>
0021
0022 #if defined(__sun) || defined(_MSC_VER)
0023 //Microsoft and solaris definition of cmath does not include math.h which has the definitions of numerical constants
0024 #include <math.h>
0025 #endif
0026
0027
0028 #ifdef HAVE_NO_EXPM1
0029 // needed to implement expm1
0030 #include <limits>
0031 #endif
0032
0033
0034 #ifndef M_PI
0035
0036 #define M_PI       3.14159265358979323846264338328      // Pi
0037 #endif
0038
0039 #ifndef M_PI_2
0040 #define M_PI_2     1.57079632679489661923132169164      // Pi/2
0041 #endif
0042
0043 #ifndef M_PI_4
0044 #define M_PI_4     0.78539816339744830961566084582      // Pi/4
0045 #endif
0046
0047 /**
0048    \namespace ROOT
0049    Namespace for new ROOT classes and functions
0050  */
0051
0052 namespace ROOT {
0053
0054 /**
0055 \namespace Math
0056 Namespace for new Math classes and functions.
0057 See the \ref Math "Math Libraries" page for a detailed description.
0058 */
0059
0060 namespace Math {
0061 // Enable Vc/VecCore template instantiations to replace std math functions.
0062 //
0063 // Vc declares `std::sqrt(Vc-type)`. To use this for Vc-`SCALAR`s, the call
0064 // to `sqrt()` must only be resolved at the template instantiation time, when
0065 // the Vc headers are guaranteed to be included, and thus its `sqrt()`
0066 // overloads have been declared.
0067 // The trick is to keep sqrt() dependent (on its argument type) by making it
0068 // an unqualified name. The `std::` of `std::sqrt()` makes it a qualified
0069 // name, so the code here has to use `sqrt()`, not `std::sqrt()`. To still
0070 // find `std::sqrt()` we pull `std::sqrt()` into the surrounding namespace.
0071 //
0072 // We don't want to use 'using namespace std' because it would pollute the including headers.
0073 using std::atan2;
0074 using std::cos;
0075 using std::cosh;
0076 using std::exp;
0077 using std::floor;
0078 using std::log;
0079 using std::pow;
0080 using std::sin;
0081 using std::sinh;
0082 using std::sqrt;
0083 using std::tan;
0084
0085 /**
0086     Mathematical constants
0087 */
0088 inline double Pi()
0089 {
0090    return M_PI;
0091    }
0092
0093    /**
0094        declarations for functions which are not implemented by some compilers
0095    */
0096
0097    /// log(1+x) with error cancelation when x is small
0098    inline double log1p(double x)
0099    {
0100 #ifndef HAVE_NO_LOG1P
0101    return ::log1p(x);
0102 #else
0103    // if log1p is not in c math library
0104   volatile double y;
0105   y = 1 + x;
0106   return std::log(y) - ((y-1)-x)/y ;  /* cancels errors with IEEE arithmetic */
0107 #endif
0108 }
0109 /// exp(x) -1 with error cancellation when x is small
0110 inline double expm1( double x) {
0111 #ifndef HAVE_NO_EXPM1
0112    return ::expm1(x);
0113 #else
0114    // compute using taylor expansion until difference is less than epsilon
0115    // use for values smaller than 0.5 (for larger (exp(x)-1 is fine
0116    if (std::abs(x) < 0.5)
0117    {
0118        // taylor series S = x + (1/2!) x^2 + (1/3!) x^3 + ...
0119
0120       double i = 1.0;
0121       double sum = x;
0122       double term = x / 1.0;
0123       do {
0124          i++ ;
0125          term *= x/i;
0126          sum += term;
0127       }
0128       while (std::abs(term) > std::abs(sum) * std::numeric_limits<double>::epsilon() ) ;
0129
0130       return sum ;
0131    }
0132    else
0133    {
0134       return std::exp(x) - 1;
0135    }
0136 #endif
0137 }
0138
0139    } // end namespace Math
0140
0141 } // end namespace ROOT
0142
0143
0144
0145
0146
0147 #endif /* ROOT_Math_Math */