// $Id$
#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 <boost/type_traits/integral_constant.hpp>
#include <limits>
#include <cmath>


namespace utl {

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


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


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


  template<typename T>
  inline
  int
  Sign(const T x)
  {
    return Sign(x, boost::integral_constant<bool, std::numeric_limits<T>::is_signed>());
  }


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

}


#include <utl/NormalDistribution.h>


namespace utl {

  /** LogGamma function
    \image html utl_LogGamma.png
    \ingroup math
  */

  double LogGamma(const double x);


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


  /**
    \brief Incomplete gamma P(a, x) function

    Taken from Numerical Recipes in C

    \ingroup math
  */

  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);
  }


  /**
    \brief Incomplete gamma Q(a, x) = 1 - P(a, x) function

    Taken from Numerical Recipes in C

    \ingroup math
  */

  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);
  }


  /** \ingroup math */

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


  /** \brief Fermi function
    \image html utl_Fermi.png

    Fermi function is defined as
    \f[
      {\rm F}(x) = \frac1{1+\exp\left(\frac{x-x_0}\sigma\right)}
    \f]

    \author Darko Veberic
    \ingroup math
  */

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


  /** Propagate axis error into theta and phi */

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


  /** Accurate solution to quadratic equation
    a x^2 + b y + c = 0
    For real solutions test return value for nan, ie !std::isnan(x.first)
  */

  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