Utilities/Math/Math.h

Go to the documentation of this file.

#ifndef _utl_Math_h_
#define _utl_Math_h_

#include <utl/MathConstants.h>
#include <utl/Vector.h>
#include <utl/Moyal.h>
#include <utl/LambertW.h>
#include <utl/Kolmogorov.h>

#include <limits>
#include <numeric>
#include <cmath>
#include <iterator>


namespace utl {

  template<typename T>
  inline
  constexpr
  T
  Sqr(const T x)
  {
    return x*x;
  }


  template<typename T>
  inline
  constexpr
  int
  Sign(const T x, const boost::integral_constant<bool, false>)
  {
    return T(0) < x;
  }


  template<typename T>
  inline
  constexpr
  int
  Sign(const T x, const boost::integral_constant<bool, true>)
  {
    return (T(0) < x) - (x < T(0));
  }


  template<typename T>
  inline
  constexpr
  int
  Sign(const T x)
  {
    // Sign(x) returns -1, 0, or 1 in case x is smaller, equal, or greater than 0, respectively
    // much like the new spaceship operator <=>
    return Sign(x, boost::integral_constant<bool, std::numeric_limits<T>::is_signed>());
  }



  inline
  double
  NormalizeAngleMinusPiPi(const double x)
  {
    return x - kTwoPi * std::floor(x/kTwoPi + 0.5);
  }



  inline
  double
  NormalizeAngleZero2Pi(const double x)
  {
    return x - kTwoPi * std::floor(x/kTwoPi);
  }


  template<typename T, std::size_t n>
  inline
  constexpr
  std::size_t
  Length(const T (&)[n])
    noexcept
  {
    return n;
  }


  inline
  constexpr
  int
  FloorDiv(const int num, const int den)
  {
    // equivalent to int(floor(double(num)/den)) but without floating point operations
    return num >= 0 ? num / den : (num - (den - 1)) / den;
  }


  namespace {

    namespace Detail {

      template<typename Container>
      inline
      double
      EvalPoly(const Container& a, const double x)
      {
        auto it = std::crbegin(a);
        const auto end = std::crend(a);
        if (it == end)
          return std::numeric_limits<double>::quiet_NaN();
        double sum = *it;
        for (++it; it != end; ++it) {
          sum *= x;
          sum += *it;
        }
        return sum;
      }

    }

  }


  template<typename T>
  inline
  double
  EvalPoly(const T& a, const double x)
  {
    return Detail::EvalPoly(a, x);
  }

  // std::initializer_list requires explicit matching ...
  template<typename T>
  inline
  double
  EvalPoly(const std::initializer_list<T>& a, const double x)
  {
    return Detail::EvalPoly(a, x);
  }

}


#include <utl/NormalDistribution.h>


namespace utl {

  inline double LogGamma(const double x) { return std::lgamma(x); }


  double IncompleteGammaPSeries(const double a, const double x);
  double IncompleteGammaPCF(const double a, const double x);


  inline
  double
  IncompleteGammaP(const double a, const double x)
  {
    // Use the series representation or the continued fractions
    return  (x < a+1) ? IncompleteGammaPSeries(a, x) : 1 - IncompleteGammaPCF(a, x);
  }


  inline
  double
  IncompleteGammaQ(const double a, const double x)
  {
    // Use the series representation or the continued fractions
    return  (x < a+1) ? 1 - IncompleteGammaPSeries(a, x) : IncompleteGammaPCF(a, x);
  }


  inline
  double
  Chi2Probability(const double chi2, const double ndof)
  {
    return IncompleteGammaQ(0.5*ndof, 0.5*chi2);
  }


  double IncompleteBeta(const double a, const double b, const double x);


  inline
  double
  Fermi(const double x, const double x0, const double sigma)
  {
    return 1/(1 + std::exp((x - x0)/sigma));
  }


  void PropagateAxisErrors(const utl::Vector::Triple& u_v_w,
                           const utl::Vector::Triple& sigma_u2_uv_v2,
                           double& thetaError, double& phiError, double& thetaPhiCorrelation);


  inline
  std::pair<double, double>
  SolveQuadraticEquation(const double a, const double b, const double c)
  {
    const double b2 = b * b;
    const double fac = 4*a*c;
    const double d = b2 - fac;
    const double sd = sqrt(d);
    if (b2 < 1e5*fac) {
      const double ta = 0.5 / a;
      return std::make_pair((-b + sd) * ta, (-b - sd) * ta);
    }
    if (b > 0) {
      const double t = -(b + sd);
      return std::make_pair(2 * c / t, t / (2 * a));
    }
    const double t = -b + sd;
    return std::make_pair(t / (2 * a), 2 * c / t);
  }

}


#endif