ATOOLS/Math/MathTools.H

0001 #ifndef ATOOLS_Math_MathTools_H
0002 #define ATOOLS_Math_MathTools_H
0003 /*  Declarations for discrete functions  */
0004
0005 #include <cmath>
0006 #include <cstdlib>
0007
0008 #include "ATOOLS/Math/MyComplex.H"
0009
0010 namespace ATOOLS {
0011
0012   template <class Type> const Type &Min(const Type &a,const Type &b)
0013   { return a<b?a:b; }
0014   template <class Type> const Type &Max(const Type &a,const Type &b)
0015   { return a>b?a:b; }
0016
0017   template <class Type> Type &Min(Type &a,Type &b)
0018   { return a<b?a:b; }
0019   template <class Type> Type &Max(Type &a,Type &b)
0020   { return a>b?a:b; }
0021
0022   inline double Accu() {return 1.e-12;}
0023   inline double SqrtAccu() {return 1.e-6;}
0024
0025   inline int    Sign(const int& a)    { return a<0?-1:1;       }
0026   inline double Sign(const double& a) { return a<0.0?-1.0:1.0; }
0027
0028   inline double Theta(const double &a) { return a<0.0?0.0:1.0; }
0029
0030   inline int    iabs(const int& a)    { return a<0?-a:a;   }
0031   inline double dabs(const double& a) { return a<0.0?-a:a; }
0032
0033   template <typename Scalar>
0034   inline Scalar sqr(const Scalar &x)   { return x*x; }
0035   template <typename Scalar> inline std::complex<Scalar>
0036   csqr(const std::complex<Scalar> &x) { return x*x; }
0037
0038
0039   inline int IsZero(const double &a,const double &crit)
0040   { return dabs(a)<crit?1:0; }
0041   inline int IsZero(const Complex &a,const double &crit)
0042   { return std::abs(a)<crit?1:0; }
0043
0044   inline int IsEqual(const double &a,const double &b)
0045   {
0046     if (a==0. && b==0.) return 1;
0047     return (dabs(a-b)/(dabs(a)+dabs(b))<Accu()) ? 1 : 0;
0048   }
0049   inline int IsEqual(const double &a,const double &b,const double &crit)
0050   {
0051     if (a==0. && b==0.) return 1;
0052     return (dabs(a-b)/(dabs(a)+dabs(b))<crit) ? 1 : 0;
0053   }
0054   inline int IsEqual(const Complex &a,const Complex &b)
0055   {
0056     if (a==Complex(0.,0.) && b==Complex(0.,0.)) return 1;
0057     return (std::abs(a-b)/(std::abs(a)+std::abs(b))<Accu()) ? 1 : 0;
0058   }
0059   inline Complex csqrt(const double &d)
0060   {
0061     if (d<0) return Complex(0.,std::sqrt(-d));
0062     return std::sqrt(d);
0063   }
0064
0065   inline double intpow(const double &a,int b)
0066   {
0067     double apb((b%2)?(b>0?a:1./a):1.);
0068     if (b>0) { while (b>1) { apb*=a*a; b-=2; } }
0069     else { while (b<-1) { apb/=a*a; b+=2; } }
0070     return apb;
0071   }
0072
0073 #define GAMMA_E 0.5772156649015328606
0074
0075   double Gammln(double xx);
0076
0077   double ReIncompleteGamma0(double x,double prec=1.e-6);
0078
0079   double  DiLog(double x);
0080   Complex DiLog(const Complex& x);
0081
0082   int Factorial(const int n);
0083
0084   double ExpIntegral(int n, double x);
0085
0086   double evaluate_polynomial (double x);
0087   long double cyl_bessel_0 (long double x);
0088   long double cyl_bessel_1 (long double x);
0089   long double cyl_bessel_2 (long double x);
0090
0091   template<typename Scalar> inline bool IsNan(const Scalar& x);
0092   template<typename Scalar> inline bool IsBad(const Scalar& x);
0093   template<typename Scalar> inline bool IsZero(const Scalar& x);
0094   template<typename Scalar> inline Scalar Abs(const Scalar& x);
0095   template<typename Scalar> inline Scalar Abs(const std::complex<Scalar>& x);
0096
0097   template<> inline bool IsNan<double>(const double& x) {
0098     return std::isnan(x)||std::isnan(-x);
0099   }
0100   template<> inline bool IsBad<double>(const double& x) {
0101     return IsNan(x)||std::isinf(x)||std::isinf(-x);
0102   }
0103   template<> inline bool IsZero<double>(const double& x) {
0104     return dabs(x)<Accu()?1:0;
0105   }
0106   template<> inline double Abs<double>(const double& x) {
0107     return x>0.0?x:-x;
0108   }
0109
0110   template<> inline bool IsNan<long double>(const long double& x) {
0111     return std::isnan(x)||std::isnan(-x);
0112   }
0113   template<> inline bool IsBad<long double>(const long double& x) {
0114     return IsNan(x)||std::isinf(x)||std::isinf(-x);
0115   }
0116   template<> inline bool IsZero<long double>(const long double& x) {
0117     return dabs(x)<Accu()?1:0;
0118   }
0119   template<> inline long double Abs<long double>(const long double& x) {
0120     return x>0.0?x:-x;
0121   }
0122
0123
0124   template<> inline bool IsNan<Complex>(const Complex& x) {
0125     return (std::isnan(real(x)) || std::isnan(imag(x)) ||
0126         std::isnan(-real(x)) || std::isnan(-imag(x)));
0127   }
0128   template<> inline bool IsBad<Complex>(const Complex& x) {
0129     return IsNan(x)||std::isinf(real(x))||std::isinf(imag(x))
0130       ||std::isinf(-real(x))||std::isinf(-imag(x));
0131   }
0132   template<> inline bool IsZero<Complex>(const Complex& x) {
0133     return std::abs(x)<Accu()?1:0;
0134   }
0135   template<> inline double Abs<double>(const Complex& x) {
0136     return std::abs(x);
0137   }
0138
0139   template<> inline bool IsNan<std::complex<long double> >
0140   (const std::complex<long double>& x) {
0141     return (std::isnan(real(x)) || std::isnan(imag(x)) ||
0142         std::isnan(-real(x)) || std::isnan(-imag(x)));
0143   }
0144   template<> inline bool IsBad<std::complex<long double> >
0145   (const std::complex<long double>& x) {
0146     return IsNan(x)||std::isinf(real(x))||std::isinf(imag(x))
0147       ||std::isinf(-real(x))||std::isinf(-imag(x));
0148   }
0149   template<> inline bool IsZero<std::complex<long double> >
0150   (const std::complex<long double>& x) {
0151     return std::abs(x)<Accu()?1:0;
0152   }
0153   template<> inline long double Abs<long double>
0154   (const std::complex<long double>& x) {
0155     return std::abs(x);
0156   }
0157
0158
0159
0160   template<class T1, class T2>
0161   struct promote_trait {
0162   };
0163
0164 #define DECLARE_PROMOTE(A,B,C)           \
0165   template<> struct promote_trait<A,B> { \
0166     typedef C T_promote;                 \
0167   }
0168
0169   DECLARE_PROMOTE(double,Complex,Complex);
0170   DECLARE_PROMOTE(Complex,double,Complex);
0171   DECLARE_PROMOTE(int,double,double);
0172   DECLARE_PROMOTE(double,int,double);
0173   DECLARE_PROMOTE(int,Complex,Complex);
0174   DECLARE_PROMOTE(Complex,int,Complex);
0175
0176   DECLARE_PROMOTE(double,double,double);
0177   DECLARE_PROMOTE(Complex,Complex,Complex);
0178
0179   DECLARE_PROMOTE(long double,std::complex<long double>,
0180           std::complex<long double>);
0181   DECLARE_PROMOTE(std::complex<long double>,long double,
0182           std::complex<long double>);
0183   DECLARE_PROMOTE(int,long double,long double);
0184   DECLARE_PROMOTE(long double,int,long double);
0185   DECLARE_PROMOTE(int,std::complex<long double>,std::complex<long double>);
0186   DECLARE_PROMOTE(std::complex<long double>,int,std::complex<long double>);
0187
0188   DECLARE_PROMOTE(long double,long double,long double);
0189   DECLARE_PROMOTE(std::complex<long double>,std::complex<long double>,
0190           std::complex<long double>);
0191
0192 #define PROMOTE(TYPE1,TYPE2) typename promote_trait<TYPE1,TYPE2>::T_promote
0193
0194   /*!
0195     \file
0196     \brief contains a collection of simple mathematical functions
0197   */
0198
0199   /*!
0200     \fn inline Type Min(Type a, Type b)
0201     \brief returns the minimum of two numbers
0202   */
0203
0204   /*!
0205     \fn inline Type Max(Type a, Type b)
0206     \brief returns the maximum of two numbers
0207   */
0208
0209   /*!
0210     \fn  inline int         Sign(const int& a) {return (a<0) ? -1 : 1;}
0211     \brief returns the sign of the argument
0212   */
0213
0214   /*!
0215     \fn  inline int         iabs(const int& a) {return a>0 ? a : -a;}
0216     \brief returns the absolute value of the argument
0217   */
0218
0219   /*!
0220     \fn  inline double      dabs(const double& a) {return a>0 ? a : -a;}
0221     \brief returns the absolute value of the argument
0222   */
0223
0224   /*!
0225     \fn  inline double      sqr(double x) {return x*x;}
0226     \brief returns the argument squared
0227   */
0228
0229   /*!
0230     \fn  inline double      Accu() {return 1.e-12;};
0231     \brief returns a (platform dependent) precission, default is \f$1^{-12}\f$
0232   */
0233
0234   /*!
0235     \fn  inline int IsZero(const double a)
0236     \brief returns \em true if argument is smaller than Accu()
0237   */
0238
0239   /*!
0240     \fn  inline int IsZero(const Complex& a)
0241     \brief  returns \em true if argument is smaller than Accu()
0242   */
0243
0244   /*!
0245     \fn  inline int IsEqual(const double a,const double b)
0246     \brief  returns \em true if arguments are equal (compared to Accu())
0247   */
0248
0249   /*!
0250     \fn  inline int IsEqual(const Complex& a,const Complex& b)
0251     \brief  returns \em true if arguments are equal (compared to Accu())
0252   */
0253
0254   /*!
0255     \fn  inline Complex csqrt(const double d)
0256     \brief returns the complex root of a (possibly negative) float or double variable
0257   */
0258
0259   /*!
0260     \fn  inline Complex csqr(Complex x)
0261     \brief returns the argument squared
0262   */
0263
0264   /*!
0265     \fn    double Gammln(double xx)
0266     \brief calculates the logarithm of the Gammafunction
0267   */
0268
0269   /*!
0270     \fn    double ReIncompleteGamma0(double xx)
0271     \brief calculates the real part of the incomplete Gammafunction.
0272   */
0273
0274   /*!
0275     \fn    double DiLog(double x)
0276     \brief calculates the real part of Li_2(x).
0277   */
0278
0279 }
0280
0281 #endif