JMath/JMathSupportkit.hh

Go to the documentation of this file.

#ifndef __JMATHSUPPORTKIT__
#define __JMATHSUPPORTKIT__

#include <limits>
#include <cmath>

#include "JTools/JConstants.hh"
#include "JLang/JException.hh"


namespace JMATH {}
namespace JPP { using namespace JMATH; }

namespace JMATH {

  using JLANG::JValueOutOfRange;

  
  /**
   * Gauss function (normalised to 1 at x = 0).
   *
   * \param  x                    x
   * \param  sigma                sigma
   * \return                      function value
   */
  inline double gauss(const double x, const double sigma)
  {
    const double u = x / sigma;

    if (fabs(u) < 10.0)
      return exp(-0.5*u*u);
    else
      return 0.0;
  }


  /**
   * Gauss function (normalised to 1 at x = x0).
   *
   * \param  x                    x
   * \param  x0                   central value
   * \param  sigma                sigma
   * \return                      function value
   */
  inline double gauss(const double x, const double x0, const double sigma)
  {
    return gauss(x - x0, sigma);
  }


  /**
   * Normalised Gauss function.
   *
   * \param  x                    x
   * \param  sigma                sigma
   * \return                      function value
   */
  inline double Gauss(const double x, const double sigma)
  {
    return gauss(x, sigma) / sqrt(2.0*JTOOLS::PI) / sigma;
  }


  /**
   * Normalised Gauss function.
   *
   * \param  x                    x
   * \param  x0                   central value
   * \param  sigma                sigma
   * \return                      function value
   */
  inline double Gauss(const double x, const double x0, const double sigma)
  {
    return Gauss(x - x0, sigma);
  }


  /**
   * Incomplete gamma function.
   *
   * \param  a                    a
   * \param  x                    x
   * \return                      function value
   */
  inline double Gamma(const double a, const double x)
  {
    using namespace std;

    const int max = 100;

    if (x <  0.0) { THROW(JValueOutOfRange, "x <  0 " << x); }
    if (a <= 0.0) { THROW(JValueOutOfRange, "a <= 0 " << a); }
 
    const double gln = lgamma(a);

    if (x <  a + 1.0) {

      if (x <= 0.0) { 
        THROW(JValueOutOfRange, "x <= 0 " << x);
      }

      double ap  = a;
      double sum = 1.0 /a;
      double del = sum;

      for (int i = 1; i != max; ++i) {

        ap  += 1.0;
        del *= x/ap;
        sum += del;

        if (fabs(del) < fabs(sum)*numeric_limits<double>::epsilon()) {
          return sum*exp(-x + a*log(x) - gln);
        }
      }

    } else {

      double b = x + 1.0 - a;
      double c = numeric_limits<double>::epsilon() / numeric_limits<double>::min();
      double d = 1.0 / b;
      double h = d;

      for (int i = 1; i != max; ++i) {

        const double an = -i * (i-a);

        b += 2.0;
        d  = an*d + b;

        if (fabs(d) < numeric_limits<double>::min()) {
          d = numeric_limits<double>::min();
        }

        c  = b + an/c;

        if (fabs(c) < numeric_limits<double>::min()) {
          c = numeric_limits<double>::min();
        }

        d  = 1.0/d;

        const double del = d*c;
        
        h *= del;
        
        if (fabs(del - 1.0) < numeric_limits<double>::epsilon()) {
          return 1.0 - exp(-x + a*log(x) - gln) * h;
        }
      }

      THROW(JValueOutOfRange, "a " << a);
    }

    return 0.0;
  }
}

#endif