ThePEG  1.8.0
Maths.h
1 // -*- C++ -*-
2 //
3 // Maths.h is a part of ThePEG - Toolkit for HEP Event Generation
4 // Copyright (C) 1999-2011 Leif Lonnblad
5 //
6 // ThePEG is licenced under version 2 of the GPL, see COPYING for details.
7 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
8 //
9 #ifndef ThePEG_Math_H
10 #define ThePEG_Math_H
11 
12 #include <cmath>
13 
14 namespace ThePEG {
15 
18 namespace Math {
19 
24 struct MathType {};
25 
27 double gamma(double);
28 
30 double lngamma(double);
31 
33 double atanh(double);
34 
37 double exp1m(double x);
38 
41 double log1m(double);
42 
44 double powi(double x, int p);
45 
47 inline double pIntegrate(double p, double xl, double xu) {
48  return p == -1.0? log(xu/xl): (pow(xu, p + 1.0) - pow(xl, p + 1.0))/(p + 1.0);
49 }
50 
52 inline double pIntegrate(int p, double xl, double xu) {
53  return p == -1? log(xu/xl): (powi(xu, p + 1) - powi(xl, p + 1))/double(p + 1);
54 }
55 
59 inline double pXIntegrate(double e, double xl, double dx) {
60  return e == 0.0? log1m(-dx/xl): -pow(xl, e)*exp1m(e*log1m(-dx/xl))/e;
61 }
62 
64 inline double pGenerate(double p, double xl, double xu, double rnd) {
65  return p == -1.0? xl*pow(xu/xl, rnd):
66  pow((1.0 - rnd)*pow(xl, p + 1.0) + rnd*pow(xu, p + 1.0), 1.0/(1.0 + p));
67 }
68 
70 inline double pGenerate(int p, double xl, double xu, double rnd) {
71  return p == -1? xl*pow(xu/xl, rnd):
72  pow((1.0 - rnd)*powi(xl, p + 1) + rnd*powi(xu, p + 1), 1.0/double(1 + p));
73 }
74 
82 inline double pXGenerate(double e, double xl, double dx, double rnd) {
83  return e == 0.0? -xl*exp1m(rnd*log1m(-dx/xl)):
84  -exp1m(log1m(rnd*exp1m(e*log1m(-dx/xl)))/e)*xl;
85 }
86 
88 template <typename FloatType>
89 inline double relativeError(FloatType x, FloatType y) {
90  return ( x == y ? 0.0 : double((x - y)/(abs(x) + abs(y))) );
91 }
92 
94 template <typename T>
95 inline T absmin(const T & x, const T & y) {
96  return abs(x) < abs(y)? x: y;
97 }
98 
100 template <typename T>
101 inline T absmax(const T & x, const T & y) {
102  return abs(x) > abs(y)? x: y;
103 }
104 
108 template <typename T, typename U>
109 inline T sign(T x, U y) {
110  return y > U()? abs(x): -abs(x);
111 }
112 
118 template <int N, bool Inv>
119 struct Power: public MathType {};
120 
124 template <int N>
125 struct Power<N,false> {
127  static double pow(double x) { return x*Power<N-1,false>::pow(x); }
128 };
129 
133 template <int N>
134 struct Power<N,true> {
136  static double pow(double x) { return Power<N+1,true>::pow(x)/x; }
137 };
138 
142 template <>
143 struct Power<0,true> {
145  static double pow(double) { return 1.0; }
146 };
147 
151 template <>
152 struct Power<0,false> {
154  static double pow(double) { return 1.0; }
155 };
157 
160 template <int N>
161 inline double Pow(double x) { return Power<N, (N < 0)>::pow(x); }
162 
166 namespace Functions {
167 
169 template <int N>
170 struct PowX: public MathType {
171 
173  static double primitive(double x) {
174  return Pow<N+1>(x)/double(N+1);
175  }
176 
178  static double integrate(double x0, double x1) {
179  return primitive(x1) - primitive(x0);
180  }
181 
184  static double generate(double x0, double x1, double R) {
185  return pow(primitive(x0) + R*integrate(x0, x1), 1.0/double(N+1));
186  }
187 
188 };
189 
195 template <>
196 inline double PowX<1>::generate(double x0, double x1, double R) {
197  return std::sqrt(x0*x0 + R*(x1*x1 - x0*x0));
198 }
199 
203 template <>
204 inline double PowX<0>::generate(double x0, double x1, double R) {
205  return x0 + R*(x1 - x0);
206 }
207 
211 template<>
212 inline double PowX<-1>::primitive(double x) {
213  return log(x);
214 }
215 
219 template <>
220 inline double PowX<-1>::integrate(double x0, double x1) {
221  return log(x1/x0);
222 }
223 
227 template <>
228 inline double PowX<-1>::generate(double x0, double x1, double R) {
229  return x0*pow(x1/x0, R);
230 }
231 
235 template <>
236 inline double PowX<-2>::generate(double x0, double x1, double R) {
237  return x0*x1/(x1 - R*(x1 - x0));
238 }
239 
243 template <>
244 inline double PowX<-3>::generate(double x0, double x1, double R) {
245  return x0*x1/std::sqrt(x1*x1 - R*(x1*x1 - x0*x0));
246 }
247 
254 template <int N>
255 struct Pow1mX: public MathType {
256 
258  static double primitive(double x) {
259  return -PowX<N>::primitive(1.0 - x);
260  }
261 
263  static double integrate(double x0, double x1) {
264  return PowX<N>::integrate(1.0 - x1, 1.0 - x0);
265  }
266 
269  static double generate(double x0, double x1, double R) {
270  return 1.0 - PowX<N>::generate(1.0 - x1, 1.0 - x0, R);
271  }
272 
273 };
274 
276 struct InvX1mX: public MathType {
277 
279  static double primitive(double x) {
280  return log(x/(1.0 - x));
281  }
282 
284  static double integrate(double x0, double x1) {
285  return log(x1*(1.0 - x0)/(x0*(1.0 - x1)));
286  }
287 
290  static double generate(double x0, double x1, double R) {
291  double r = pow(x1*(1.0 - x0)/(x0*(1.0 - x1)), R)*x0/(1.0 - x0);
292  return r/(1.0 + r);
293  }
294 
295 };
296 
298 struct ExpX: public MathType {
299 
301  static double primitive(double x) {
302  return exp(x);
303  }
304 
306  static double integrate(double x0, double x1) {
307  return exp(x1) - exp(x0);
308  }
309 
312  static double generate(double x0, double x1, double R) {
313  return log(exp(x0) + R*(exp(x1) - exp(x0)));
314  }
315 
316 };
317 
320 template <int N, int D>
321 struct FracPowX: public MathType {
322 
324  static double primitive(double x) {
325  double r = double(N)/double(D) + 1.0;
326  return pow(x, r)/r;
327  }
328 
330  static double integrate(double x0, double x1) {
331  return primitive(x1) - primitive(x0);
332  }
333 
336  static double generate(double x0, double x1, double R) {
337  double r = double(N)/double(D) + 1.0;
338  return pow(primitive(x0) + R*integrate(x0, x1), 1.0/r);
339  }
340 
341 };
342 
343 }
344 
345 }
346 
347 }
348 
349 #endif /* ThePEG_Math_H */
Class corresponding to functions of the form with integer N.
Definition: Maths.h:170
static double generate(double x0, double x1, double R)
Sample a distribution in a given interval given a flat random number R in the interval ]0...
Definition: Maths.h:290
Class corresponding to functions of the form with integer N.
Definition: Maths.h:255
static double integrate(double x0, double x1)
Integrate function in a given interval.
Definition: Maths.h:306
static double integrate(double x0, double x1)
Integrate function in a given interval.
Definition: Maths.h:263
static double primitive(double x)
The primitive function.
Definition: Maths.h:279
Class corresponding to functions of the form .
Definition: Maths.h:276
double pIntegrate(double p, double xl, double xu)
Return the integral of between xl and xu.
Definition: Maths.h:47
static double pow(double x)
Member for the power.
Definition: Maths.h:127
This is the main namespace within which all identifiers in ThePEG are declared.
Definition: FactoryBase.h:28
static double pow(double)
Member for the power.
Definition: Maths.h:145
double exp1m(double x)
Return , with highest possible precision for .
static double pow(double x)
Member for the power.
Definition: Maths.h:136
T absmin(const T &x, const T &y)
Return x if |x|<|y|, else return y.
Definition: Maths.h:95
double pXIntegrate(double e, double xl, double dx)
Return the integral of between xl and xl+dx with highest possible precision for and/or ...
Definition: Maths.h:59
double log1m(double)
Return , with highest possible precision for .
Templated class for calculating integer powers.
Definition: Maths.h:119
static double integrate(double x0, double x1)
Integrate function in a given interval.
Definition: Maths.h:284
static double primitive(double x)
The primitive function.
Definition: Maths.h:173
Class corresponding to functions of the form with integer N and D.
Definition: Maths.h:321
static double primitive(double x)
The primitive function.
Definition: Maths.h:301
double relativeError(FloatType x, FloatType y)
Returns (x - y)/(|x| + |y|).
Definition: Maths.h:89
static double integrate(double x0, double x1)
Integrate function in a given interval.
Definition: Maths.h:330
double powi(double x, int p)
Return x rased to the integer power p, using recursion.
double gamma(double)
The gamma function.
double pXGenerate(double e, double xl, double dx, double rnd)
Generate an x between xl and xl + dx distributed as with highest possible precision for and/or * ...
Definition: Maths.h:82
static double pow(double)
Member for the power.
Definition: Maths.h:154
double pGenerate(double p, double xl, double xu, double rnd)
Generate an x between xl and xu distributed as .
Definition: Maths.h:64
Class corresponding to functions of the form .
Definition: Maths.h:298
static double primitive(double x)
The primitive function.
Definition: Maths.h:258
static double generate(double x0, double x1, double R)
Sample a distribution in a given interval given a flat random number R in the interval ]0...
Definition: Maths.h:269
MathType is an empty non-polymorphic base class for all mathematical function types.
Definition: Maths.h:24
static double primitive(double x)
The primitive function.
Definition: Maths.h:324
static double integrate(double x0, double x1)
Integrate function in a given interval.
Definition: Maths.h:178
double atanh(double)
Return .
static double generate(double x0, double x1, double R)
Sample a distribution in a given interval given a flat random number R in the interval ]0...
Definition: Maths.h:312
T sign(T x, U y)
Transfer the sign of the second argument to the first.
Definition: Maths.h:109
static double generate(double x0, double x1, double R)
Sample a distribution in a given interval given a flat random number R in the interval ]0...
Definition: Maths.h:336
double lngamma(double)
The log of the gamma function.
double Pow(double x)
Templated function to calculate integer powers known at compile-time.
Definition: Maths.h:161
static double generate(double x0, double x1, double R)
Sample a distribution in a given interval given a flat random number R in the interval ]0...
Definition: Maths.h:184
T absmax(const T &x, const T &y)
Return x if |x|>|y|, else return y.
Definition: Maths.h:101