ThePEG/Utilities/Maths.h

0001 // -*- C++ -*-
0002 //
0003 // Maths.h is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 1999-2019 Leif Lonnblad
0005 //
0006 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef ThePEG_Math_H
0010 #define ThePEG_Math_H
0011
0012 #include <cmath>
0013
0014 namespace ThePEG {
0015
0016 /** The Math namespace includes the declaration of some useful
0017  *  mathematical functions. */
0018 namespace Math {
0019
0020 /**
0021  * MathType is an empty non-polymorphic base class for all
0022  * mathematical function types.
0023  */
0024 struct MathType {};
0025
0026 /** Return \f$1-e^x\f$, with highest possible precision for
0027  *  \f$x\rightarrow 0\f$. */
0028 double exp1m(double x);
0029
0030 /** Return \f$1\log(1-x)\f$, with highest possible precision for
0031  *  \f$x\rightarrow 0\f$. */
0032 double log1m(double);
0033
0034 /** Return x rased to the integer power p, using recursion. */
0035 double powi(double x, int p);
0036
0037 /** Return the integral of \f$x^p dx\f$ between xl and xu. */
0038 inline double pIntegrate(double p, double xl, double xu) {
0039   return p == -1.0? log(xu/xl): (pow(xu, p + 1.0) - pow(xl, p + 1.0))/(p + 1.0);
0040 }
0041
0042 /** Return the integral of \f$x^p dx\f$ between xl and xu. */
0043 inline double pIntegrate(int p, double xl, double xu) {
0044   return p == -1? log(xu/xl): (powi(xu, p + 1) - powi(xl, p + 1))/double(p + 1);
0045 }
0046
0047 /** Return the integral of \f$x^{e-1} dx\f$ between xl and xl+dx with
0048  *  highest possible precision for \f$dx\rightarrow 0\f$ and/or
0049  *  \f$e\rightarrow 0\f$. */
0050 inline double pXIntegrate(double e, double xl, double dx) {
0051   return e == 0.0? log1m(-dx/xl): -pow(xl, e)*exp1m(e*log1m(-dx/xl))/e;
0052 }
0053
0054 /** Generate an x between xl and xu distributed as \f$x^p\f$. */
0055 inline double pGenerate(double p, double xl, double xu, double rnd) {
0056   return p == -1.0? xl*pow(xu/xl, rnd):
0057     pow((1.0 - rnd)*pow(xl, p + 1.0) + rnd*pow(xu, p + 1.0), 1.0/(1.0 + p));
0058 }
0059
0060 /** Generate an x between xl and xu distributed as \f$x^p\f$. */
0061 inline double pGenerate(int p, double xl, double xu, double rnd) {
0062   return p == -1? xl*pow(xu/xl, rnd):
0063     pow((1.0 - rnd)*powi(xl, p + 1) + rnd*powi(xu, p + 1), 1.0/double(1 + p));
0064 }
0065
0066 /** Generate an x between xl and xl + dx distributed as \f$x^{e-1}\f$
0067  *  with highest possible precision for\f$dx\rightarrow 0\f$ and/or *
0068  *  \f$e\rightarrow 0\f$.
0069  * @param e the parameter defining the power in \f$x^{e-1}\f$.
0070  * @param xl the lower bound of the generation interval.
0071  * @param dx the interval.
0072  * @param rnd a flat random number in the interval ]0,1[. */
0073 inline double pXGenerate(double e, double xl, double dx, double rnd) {
0074   return e == 0.0? -xl*exp1m(rnd*log1m(-dx/xl)):
0075     -exp1m(log1m(rnd*exp1m(e*log1m(-dx/xl)))/e)*xl;
0076 }
0077
0078 /** Returns (x - y)/(|x| + |y|). */
0079 template <typename FloatType>
0080 inline double relativeError(FloatType x, FloatType y) {
0081   return ( x == y ? 0.0 : double((x - y)/(abs(x) + abs(y))) );
0082 }
0083
0084 /** Return x if |x|<|y|, else return y. */
0085 template <typename T>
0086 inline T absmin(const T & x, const T & y) {
0087   return abs(x) < abs(y)? x: y;
0088 }
0089
0090 /** Return x if |x|>|y|, else return y. */
0091 template <typename T>
0092 inline T absmax(const T & x, const T & y) {
0093   return abs(x) > abs(y)? x: y;
0094 }
0095
0096 /** Transfer the sign of the second argument to the first.
0097  * @return \f$|x|\f$ if \f$y>0\f$ otherwise return \f$-|x|\f$.
0098  */
0099 template <typename T, typename U>
0100 inline T sign(T x, U y) {
0101   return y > U()? abs(x): -abs(x);
0102 }
0103
0104 /** Templated class for calculating integer powers. */
0105 //@{
0106 /**
0107  *  Struct for powers
0108  */
0109 template <int N, bool Inv>
0110 struct Power: public MathType {};
0111
0112 /**
0113  *  Struct for powers
0114  */
0115 template <int N>
0116 struct Power<N,false> {
0117   /** Member for the power*/
0118   static double pow(double x) { return x*Power<N-1,false>::pow(x); }
0119 };
0120
0121 /**
0122  *  Struct for powers
0123  */
0124 template <int N>
0125 struct Power<N,true> {
0126   /** Member for the power*/
0127   static double pow(double x) { return Power<N+1,true>::pow(x)/x; }
0128 };
0129
0130 /**
0131  *  Struct for powers
0132  */
0133 template <>
0134 struct Power<0,true> {
0135   /** Member for the power*/
0136   static double pow(double) { return 1.0; }
0137 };
0138
0139 /**
0140  *  Struct for powers
0141  */
0142 template <>
0143 struct Power<0,false> {
0144   /** Member for the power*/
0145   static double pow(double) { return 1.0; }
0146 };
0147 //@}
0148
0149 /** Templated function to calculate integer powers known at
0150  *  compile-time. */
0151 template <int N>
0152 inline double Pow(double x) { return Power<N, (N < 0)>::pow(x); }
0153
0154 /** This namespace introduces some useful function classes with known
0155  *  primitive and inverse primitive functions. Useful to sample
0156  *  corresponding distributions.*/
0157 namespace Functions {
0158
0159 /** Class corresponding to functions of the form \f$x^N\f$ with integer N. */
0160 template <int N>
0161 struct PowX: public MathType {
0162
0163   /** The primitive function. */
0164   static double primitive(double x) {
0165     return Pow<N+1>(x)/double(N+1);
0166   }
0167
0168   /** Integrate function in a given interval. */
0169   static double integrate(double x0, double x1) {
0170     return primitive(x1) - primitive(x0);
0171   }
0172
0173   /** Sample a distribution in a given interval given a flat random
0174    *  number R in the interval ]0,1[. */
0175   static double generate(double x0, double x1, double R) {
0176     return pow(primitive(x0) + R*integrate(x0, x1), 1.0/double(N+1));
0177   }
0178
0179 };
0180
0181 /** @cond TRAITSPECIALIZATIONS */
0182
0183 /**
0184  *  Template for generating according to a specific power
0185  */
0186 template <>
0187 inline double PowX<1>::generate(double x0, double x1, double R) {
0188   return std::sqrt(x0*x0 + R*(x1*x1 - x0*x0));
0189 }
0190
0191 /**
0192  *  Template for generating according to a specific power
0193  */
0194 template <>
0195 inline double PowX<0>::generate(double x0, double x1, double R) {
0196   return x0 + R*(x1 - x0);
0197 }
0198
0199 /**
0200  *  Template for generating according to a specific power
0201  */
0202 template<>
0203 inline double PowX<-1>::primitive(double x) {
0204   return log(x);
0205 }
0206
0207 /**
0208  *  Template for generating according to a specific power
0209  */
0210 template <>
0211 inline double PowX<-1>::integrate(double x0, double x1) {
0212   return log(x1/x0);
0213 }
0214
0215 /**
0216  *  Template for generating according to a specific power
0217  */
0218 template <>
0219 inline double PowX<-1>::generate(double x0, double x1, double R) {
0220   return x0*pow(x1/x0, R);
0221 }
0222
0223 /**
0224  *  Template for generating according to a specific power
0225  */
0226 template <>
0227 inline double PowX<-2>::generate(double x0, double x1, double R) {
0228   return x0*x1/(x1 - R*(x1 - x0));
0229 }
0230
0231 /**
0232  *  Template for generating according to a specific power
0233  */
0234 template <>
0235 inline double PowX<-3>::generate(double x0, double x1, double R) {
0236   return x0*x1/std::sqrt(x1*x1 - R*(x1*x1 - x0*x0));
0237 }
0238
0239 /** @endcond */
0240
0241
0242
0243 /** Class corresponding to functions of the form \f$(1-x)^N\f$
0244  *  with integer N. */
0245 template <int N>
0246 struct Pow1mX: public MathType {
0247
0248   /** The primitive function. */
0249   static double primitive(double x) {
0250     return -PowX<N>::primitive(1.0 - x);
0251   }
0252
0253   /** Integrate function in a given interval. */
0254   static double integrate(double x0, double x1) {
0255     return PowX<N>::integrate(1.0 - x1, 1.0 - x0);
0256   }
0257
0258   /** Sample a distribution in a given interval given a flat random
0259    *  number R in the interval ]0,1[. */
0260   static double generate(double x0, double x1, double R) {
0261     return 1.0 - PowX<N>::generate(1.0 - x1, 1.0 - x0, R);
0262   }
0263
0264 };
0265
0266 /** Class corresponding to functions of the form \f$1/(x(1-x))\f$ */
0267 struct InvX1mX: public MathType {
0268
0269   /** The primitive function. */
0270   static double primitive(double x) {
0271     return log(x/(1.0 - x));
0272   }
0273
0274   /** Integrate function in a given interval. */
0275   static double integrate(double x0, double x1) {
0276     return log(x1*(1.0 - x0)/(x0*(1.0 - x1)));
0277   }
0278
0279   /** Sample a distribution in a given interval given a flat random
0280    *  number R in the interval ]0,1[. */
0281   static double generate(double x0, double x1, double R) {
0282     double r = pow(x1*(1.0 - x0)/(x0*(1.0 - x1)), R)*x0/(1.0 - x0);
0283     return r/(1.0 + r);
0284   }
0285
0286 };
0287
0288 /** Class corresponding to functions of the form \f$e^x\f$ */
0289 struct ExpX: public MathType {
0290
0291   /** The primitive function. */
0292   static double primitive(double x) {
0293     return exp(x);
0294   }
0295
0296   /** Integrate function in a given interval. */
0297   static double integrate(double x0, double x1) {
0298     return exp(x1) - exp(x0);
0299   }
0300
0301   /** Sample a distribution in a given interval given a flat random
0302    *  number R in the interval ]0,1[. */
0303   static double generate(double x0, double x1, double R) {
0304     return log(exp(x0) + R*(exp(x1) - exp(x0)));
0305   }
0306
0307 };
0308
0309 /** Class corresponding to functions of the form \f$x^{N/D}\f$
0310  *  with integer N and D. */
0311 template <int N, int D>
0312 struct FracPowX: public MathType {
0313
0314   /** The primitive function. */
0315   static double primitive(double x) {
0316     double r = double(N)/double(D) + 1.0;
0317     return pow(x, r)/r;
0318   }
0319
0320   /** Integrate function in a given interval. */
0321   static double integrate(double x0, double x1) {
0322     return primitive(x1) - primitive(x0);
0323   }
0324
0325   /** Sample a distribution in a given interval given a flat random
0326    *  number R in the interval ]0,1[. */
0327   static double generate(double x0, double x1, double R) {
0328     double r = double(N)/double(D) + 1.0;
0329     return pow(primitive(x0) + R*integrate(x0, x1), 1.0/r);
0330   }
0331
0332 };
0333
0334 }
0335
0336 }
0337
0338 }
0339
0340 #endif /* ThePEG_Math_H */