00001
00002
00003
00004
00005
00006
00007
00008
00009 #ifndef ThePEG_Math_H
00010 #define ThePEG_Math_H
00011
00012 #include <cmath>
00013
00014 namespace ThePEG {
00015
00018 namespace Math {
00019
00024 struct MathType {};
00025
00027 double gamma(double);
00028
00030 double lngamma(double);
00031
00033 double atanh(double);
00034
00037 double exp1m(double x);
00038
00041 double log1m(double);
00042
00044 double powi(double x, int p);
00045
00047 inline double pIntegrate(double p, double xl, double xu) {
00048 return p == -1.0? log(xu/xl): (pow(xu, p + 1.0) - pow(xl, p + 1.0))/(p + 1.0);
00049 }
00050
00052 inline double pIntegrate(int p, double xl, double xu) {
00053 return p == -1? log(xu/xl): (powi(xu, p + 1) - powi(xl, p + 1))/double(p + 1);
00054 }
00055
00059 inline double pXIntegrate(double e, double xl, double dx) {
00060 return e == 0.0? log1m(-dx/xl): -pow(xl, e)*exp1m(e*log1m(-dx/xl))/e;
00061 }
00062
00064 inline double pGenerate(double p, double xl, double xu, double rnd) {
00065 return p == -1.0? xl*pow(xu/xl, rnd):
00066 pow((1.0 - rnd)*pow(xl, p + 1.0) + rnd*pow(xu, p + 1.0), 1.0/(1.0 + p));
00067 }
00068
00070 inline double pGenerate(int p, double xl, double xu, double rnd) {
00071 return p == -1? xl*pow(xu/xl, rnd):
00072 pow((1.0 - rnd)*powi(xl, p + 1) + rnd*powi(xu, p + 1), 1.0/double(1 + p));
00073 }
00074
00082 inline double pXGenerate(double e, double xl, double dx, double rnd) {
00083 return e == 0.0? -xl*exp1m(rnd*log1m(-dx/xl)):
00084 -exp1m(log1m(rnd*exp1m(e*log1m(-dx/xl)))/e)*xl;
00085 }
00086
00088 template <typename FloatType>
00089 inline double relativeError(FloatType x, FloatType y) {
00090 return ( x == y ? 0.0 : double((x - y)/(abs(x) + abs(y))) );
00091 }
00092
00094 template <typename T>
00095 inline T absmin(const T & x, const T & y) {
00096 return abs(x) < abs(y)? x: y;
00097 }
00098
00100 template <typename T>
00101 inline T absmax(const T & x, const T & y) {
00102 return abs(x) > abs(y)? x: y;
00103 }
00104
00108 template <typename T, typename U>
00109 inline T sign(T x, U y) {
00110 return y > U()? abs(x): -abs(x);
00111 }
00112
00118 template <int N, bool Inv>
00119 struct Power: public MathType {};
00120
00124 template <int N>
00125 struct Power<N,false> {
00127 static double pow(double x) { return x*Power<N-1,false>::pow(x); }
00128 };
00129
00133 template <int N>
00134 struct Power<N,true> {
00136 static double pow(double x) { return Power<N+1,true>::pow(x)/x; }
00137 };
00138
00142 template <>
00143 struct Power<0,true> {
00145 static double pow(double) { return 1.0; }
00146 };
00147
00151 template <>
00152 struct Power<0,false> {
00154 static double pow(double) { return 1.0; }
00155 };
00157
00160 template <int N>
00161 inline double Pow(double x) { return Power<N, (N < 0)>::pow(x); }
00162
00166 namespace Functions {
00167
00169 template <int N>
00170 struct PowX: public MathType {
00171
00173 static double primitive(double x) {
00174 return Pow<N+1>(x)/double(N+1);
00175 }
00176
00178 static double integrate(double x0, double x1) {
00179 return primitive(x1) - primitive(x0);
00180 }
00181
00184 static double generate(double x0, double x1, double R) {
00185 return pow(primitive(x0) + R*integrate(x0, x1), 1.0/double(N+1));
00186 }
00187
00188 };
00189
00195 template <>
00196 inline double PowX<1>::generate(double x0, double x1, double R) {
00197 return std::sqrt(x0*x0 + R*(x1*x1 - x0*x0));
00198 }
00199
00203 template <>
00204 inline double PowX<0>::generate(double x0, double x1, double R) {
00205 return x0 + R*(x1 - x0);
00206 }
00207
00211 template<>
00212 inline double PowX<-1>::primitive(double x) {
00213 return log(x);
00214 }
00215
00219 template <>
00220 inline double PowX<-1>::integrate(double x0, double x1) {
00221 return log(x1/x0);
00222 }
00223
00227 template <>
00228 inline double PowX<-1>::generate(double x0, double x1, double R) {
00229 return x0*pow(x1/x0, R);
00230 }
00231
00235 template <>
00236 inline double PowX<-2>::generate(double x0, double x1, double R) {
00237 return x0*x1/(x1 - R*(x1 - x0));
00238 }
00239
00243 template <>
00244 inline double PowX<-3>::generate(double x0, double x1, double R) {
00245 return x0*x1/std::sqrt(x1*x1 - R*(x1*x1 - x0*x0));
00246 }
00247
00254 template <int N>
00255 struct Pow1mX: public MathType {
00256
00258 static double primitive(double x) {
00259 return -PowX<N>::primitive(1.0 - x);
00260 }
00261
00263 static double integrate(double x0, double x1) {
00264 return PowX<N>::integrate(1.0 - x1, 1.0 - x0);
00265 }
00266
00269 static double generate(double x0, double x1, double R) {
00270 return 1.0 - PowX<N>::generate(1.0 - x1, 1.0 - x0, R);
00271 }
00272
00273 };
00274
00276 struct InvX1mX: public MathType {
00277
00279 static double primitive(double x) {
00280 return log(x/(1.0 - x));
00281 }
00282
00284 static double integrate(double x0, double x1) {
00285 return log(x1*(1.0 - x0)/(x0*(1.0 - x1)));
00286 }
00287
00290 static double generate(double x0, double x1, double R) {
00291 double r = pow(x1*(1.0 - x0)/(x0*(1.0 - x1)), R)*x0/(1.0 - x0);
00292 return r/(1.0 + r);
00293 }
00294
00295 };
00296
00298 struct ExpX: public MathType {
00299
00301 static double primitive(double x) {
00302 return exp(x);
00303 }
00304
00306 static double integrate(double x0, double x1) {
00307 return exp(x1) - exp(x0);
00308 }
00309
00312 static double generate(double x0, double x1, double R) {
00313 return log(exp(x0) + R*(exp(x1) - exp(x0)));
00314 }
00315
00316 };
00317
00320 template <int N, int D>
00321 struct FracPowX: public MathType {
00322
00324 static double primitive(double x) {
00325 double r = double(N)/double(D) + 1.0;
00326 return pow(x, r)/r;
00327 }
00328
00330 static double integrate(double x0, double x1) {
00331 return primitive(x1) - primitive(x0);
00332 }
00333
00336 static double generate(double x0, double x1, double R) {
00337 double r = double(N)/double(D) + 1.0;
00338 return pow(primitive(x0) + R*integrate(x0, x1), 1.0/r);
00339 }
00340
00341 };
00342
00343 }
00344
00345 }
00346
00347 }
00348
00349 #endif