Math.h

Go to the documentation of this file.

// @(#)root/mathcore:$Id$
// Author: L. Moneta Tue Nov 14 15:44:38 2006
 
/**********************************************************************
 *                                                                    *
 * Copyright (c) 2006  LCG ROOT Math Team, CERN/PH-SFT                *
 *                                                                    *
 *                                                                    *
 **********************************************************************/
 
// mathematical constants like Pi
 
#ifndef ROOT_Math_Math
#define ROOT_Math_Math
 
#ifdef _MSC_VER
#define _USE_MATH_DEFINES
#endif
 
#include <cmath>
 
#if defined(__sun) || defined(_MSC_VER)
//Microsoft and solaris definition of cmath does not include math.h which has the definitions of numerical constants
#include <math.h> // for M_PI
// TODO replace with std::numbers::pi once minimum version is C++20, and remove this code block
#endif
 
 
#ifdef HAVE_NO_EXPM1
// needed to implement expm1
#include <limits>
#endif
 
 
#ifndef M_PI
 
#define M_PI       3.14159265358979323846264338328      // Pi
#endif
 
#ifndef M_PI_2
#define M_PI_2     1.57079632679489661923132169164      // Pi/2
#endif
 
#ifndef M_PI_4
#define M_PI_4     0.78539816339744830961566084582      // Pi/4
#endif
 
namespace ROOT {
 
/**
\namespace Math
Namespace for new Math classes and functions.
See the \ref Math "Math Libraries" page for a detailed description.
*/
 
namespace Math {
 
/**
    Mathematical constants
*/
inline double Pi()
{
   return M_PI;
   }
 
   /**
       declarations for functions which are not implemented by some compilers
   */
 
   /// log(1+x) with error cancelation when x is small
   inline double log1p(double x)
   {
#ifndef HAVE_NO_LOG1P
   return ::log1p(x);
#else
   // if log1p is not in c math library
  volatile double y;
  y = 1 + x;
  return std::log(y) - ((y-1)-x)/y ;  /* cancels errors with IEEE arithmetic */
#endif
}
/// exp(x) -1 with error cancellation when x is small
inline double expm1( double x) {
#ifndef HAVE_NO_EXPM1
   return ::expm1(x);
#else
   // compute using taylor expansion until difference is less than epsilon
   // use for values smaller than 0.5 (for larger (exp(x)-1 is fine
   if (std::abs(x) < 0.5)
   {
       // taylor series S = x + (1/2!) x^2 + (1/3!) x^3 + ...
 
      double i = 1.0;
      double sum = x;
      double term = x / 1.0;
      do {
         i++ ;
         term *= x/i;
         sum += term;
      }
      while (std::abs(term) > std::abs(sum) * std::numeric_limits<double>::epsilon() ) ;
 
      return sum ;
   }
   else
   {
      return std::exp(x) - 1;
   }
#endif
}
 
} // end namespace Math
 
} // end namespace ROOT
 
 
 
 
 
#endif /* ROOT_Math_Math */