Maths.h

// -*- C++ -*-
//
// Maths.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ThePEG_Math_H
#define ThePEG_Math_H

#include <cmath>

namespace ThePEG {

namespace Math {

struct MathType {};

double exp1m(double x);

double log1m(double);

double powi(double x, int p);

inline double pIntegrate(double p, double xl, double xu) {
  return p == -1.0? log(xu/xl): (pow(xu, p + 1.0) - pow(xl, p + 1.0))/(p + 1.0);
}

inline double pIntegrate(int p, double xl, double xu) {
  return p == -1? log(xu/xl): (powi(xu, p + 1) - powi(xl, p + 1))/double(p + 1);
}

inline double pXIntegrate(double e, double xl, double dx) {
  return e == 0.0? log1m(-dx/xl): -pow(xl, e)*exp1m(e*log1m(-dx/xl))/e;
}

inline double pGenerate(double p, double xl, double xu, double rnd) {
  return p == -1.0? xl*pow(xu/xl, rnd):
    pow((1.0 - rnd)*pow(xl, p + 1.0) + rnd*pow(xu, p + 1.0), 1.0/(1.0 + p));
}

inline double pGenerate(int p, double xl, double xu, double rnd) {
  return p == -1? xl*pow(xu/xl, rnd):
    pow((1.0 - rnd)*powi(xl, p + 1) + rnd*powi(xu, p + 1), 1.0/double(1 + p));
}

inline double pXGenerate(double e, double xl, double dx, double rnd) {
  return e == 0.0? -xl*exp1m(rnd*log1m(-dx/xl)):
    -exp1m(log1m(rnd*exp1m(e*log1m(-dx/xl)))/e)*xl;
}

template <typename FloatType>
inline double relativeError(FloatType x, FloatType y) {
  return ( x == y ? 0.0 : double((x - y)/(abs(x) + abs(y))) );
}

template <typename T>
inline T absmin(const T & x, const T & y) {
  return abs(x) < abs(y)? x: y;
}

template <typename T>
inline T absmax(const T & x, const T & y) {
  return abs(x) > abs(y)? x: y;
}

template <typename T, typename U>
inline T sign(T x, U y) {
  return y > U()? abs(x): -abs(x);
}

template <int N, bool Inv>
struct Power: public MathType {};

template <int N>
struct Power<N,false> {
  static double pow(double x) { return x*Power<N-1,false>::pow(x); }
};

template <int N>
struct Power<N,true> {
  static double pow(double x) { return Power<N+1,true>::pow(x)/x; }
};

template <>
struct Power<0,true> {
  static double pow(double) { return 1.0; }
};

template <>
struct Power<0,false> {
  static double pow(double) { return 1.0; }
};

template <int N>
inline double Pow(double x) { return Power<N, (N < 0)>::pow(x); }

namespace Functions {

template <int N>
struct PowX: public MathType {

  static double primitive(double x) { 
    return Pow<N+1>(x)/double(N+1); 
  }

  static double integrate(double x0, double x1) {
    return primitive(x1) - primitive(x0);
  }

  static double generate(double x0, double x1, double R) {
    return pow(primitive(x0) + R*integrate(x0, x1), 1.0/double(N+1));
  }

};

template <>
inline double PowX<1>::generate(double x0, double x1, double R) {
  return std::sqrt(x0*x0 + R*(x1*x1 - x0*x0));
}

template <>
inline double PowX<0>::generate(double x0, double x1, double R) {
  return x0 + R*(x1 - x0);
}

template<>
inline double PowX<-1>::primitive(double x) {
  return log(x);
}

template <>
inline double PowX<-1>::integrate(double x0, double x1) {
  return log(x1/x0);
}

template <>
inline double PowX<-1>::generate(double x0, double x1, double R) {
  return x0*pow(x1/x0, R);
}

template <>
inline double PowX<-2>::generate(double x0, double x1, double R) {
  return x0*x1/(x1 - R*(x1 - x0));
}

template <>
inline double PowX<-3>::generate(double x0, double x1, double R) {
  return x0*x1/std::sqrt(x1*x1 - R*(x1*x1 - x0*x0));
}

template <int N>
struct Pow1mX: public MathType {

  static double primitive(double x) {
    return -PowX<N>::primitive(1.0 - x);
  }

  static double integrate(double x0, double x1) {
    return PowX<N>::integrate(1.0 - x1, 1.0 - x0);
  }

  static double generate(double x0, double x1, double R) {
    return 1.0 - PowX<N>::generate(1.0 - x1, 1.0 - x0, R);
  }

};

struct InvX1mX: public MathType {

  static double primitive(double x) {
    return log(x/(1.0 - x));
  }

  static double integrate(double x0, double x1) {
    return log(x1*(1.0 - x0)/(x0*(1.0 - x1)));
  }

  static double generate(double x0, double x1, double R) {
    double r = pow(x1*(1.0 - x0)/(x0*(1.0 - x1)), R)*x0/(1.0 - x0);
    return r/(1.0 + r);
  }

};

struct ExpX: public MathType {

  static double primitive(double x) { 
    return exp(x);
  }

  static double integrate(double x0, double x1) {
    return exp(x1) - exp(x0);
  }

  static double generate(double x0, double x1, double R) {
    return log(exp(x0) + R*(exp(x1) - exp(x0)));
  }

};  

template <int N, int D>
struct FracPowX: public MathType {

  static double primitive(double x) {
    double r = double(N)/double(D) + 1.0;
    return pow(x, r)/r;
  }

  static double integrate(double x0, double x1) {
    return primitive(x1) - primitive(x0);
  }

  static double generate(double x0, double x1, double R) {
    double r = double(N)/double(D) + 1.0;
    return pow(primitive(x0) + R*integrate(x0, x1), 1.0/r);
  }

};

}

}

}

#endif /* ThePEG_Math_H */