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 */