Gamma.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
//
// ------------------------------------------------------------
// GEANT 4 class implementation
// ------------------------------------------------------------
 
#include <cmath>
#include <string.h>
#include "Gamma.hh"
 
MyGamma::MyGamma(){}
 
MyGamma::~MyGamma(){}
 
//____________________________________________________________________________
double MyGamma::Gamma(double z)
{
  if (z <= 0)
      return 0;
 
  return std::tgamma(z);
}
 
//____________________________________________________________________________
double MyGamma::Gamma(double a, double x)
{
  // Computation of the incomplete gamma function P(a,x)
  //
  // The algorithm is based on the formulas and code as denoted in
  // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
  //
  //--- Nve 14-nov-1998 UU-SAP Utrecht
 
  if (a <= 0 || x <= 0) return 0;
 
  if (x < (a + 1))
    return GamSer(a, x);
  else
    return GamCf(a, x);
}
 
//____________________________________________________________________________
double MyGamma::GamCf(double a, double x)
{
  // Computation of the incomplete gamma function P(a,x)
  // via its continued fraction representation.
  //
  // The algorithm is based on the formulas and code as denoted in
  // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
  //
  //--- Nve 14-nov-1998 UU-SAP Utrecht
 
  int itmax = 100;  // Maximum number of iterations
  double eps = 3.e-7;  // Relative accuracy
  double fpmin = 1.e-30;  // Smallest double value allowed here
 
  if (a <= 0 || x <= 0) return 0;
 
  double gln = LnGamma(a);
  double b = x + 1 - a;
  double c = 1 / fpmin;
  double d = 1 / b;
  double h = d;
  double an, del;
  for (int i = 1; i <= itmax; i++) {
    an = double(-i) * (double(i) - a);
    b += 2;
    d = an * d + b;
    if (Abs(d) < fpmin) d = fpmin;
    c = b + an / c;
    if (Abs(c) < fpmin) c = fpmin;
    d = 1 / d;
    del = d * c;
    h = h * del;
    if (Abs(del - 1) < eps) break;
    // if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl;
  }
  double v = Exp(-x + a * Log(x) - gln) * h;
  return (1 - v);
}
 
//____________________________________________________________________________
double MyGamma::GamSer(double a, double x)
{
  // Computation of the incomplete gamma function P(a,x)
  // via its series representation.
  //
  // The algorithm is based on the formulas and code as denoted in
  // Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
  //
  //--- Nve 14-nov-1998 UU-SAP Utrecht
 
  int itmax = 100;  // Maximum number of iterations
  double eps = 3.e-7;  // Relative accuracy
 
  if (a <= 0 || x <= 0) return 0;
 
  double gln = LnGamma(a);
  double ap = a;
  double sum = 1 / a;
  double del = sum;
  for (int n = 1; n <= itmax; n++) {
    ap += 1;
    del = del * x / ap;
    sum += del;
    if (MyGamma::Abs(del) < Abs(sum * eps)) break;
    // if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl;
  }
  double v = sum * Exp(-x + a * Log(x) - gln);
  return v;
}
 
double MyGamma::LnGamma(double z)
{
  if (z <= 0)
      return 0;
 
  return std::lgamma(z);
}