// @(#)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 */