NumRecipes.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------
// File and Version Information:
//  $Id: NumRecipes.cxx,v 1.1.1.1 2005/04/21 01:17:17 zhangy Exp $
//
// Description:
//
// Environment:
//  Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Author List:
//  Bob Jacobsen, Ed Iskander
//
// Copyright Information:
//  Copyright (C) 1996  Lawrence Berkeley Laboratory
//
//------------------------------------------------------------------------
//#include "BaBar/BaBar.hh"
 
//-----------------------
// This Class's Header --
//-----------------------
#include "ProbTools/NumRecipes.h"
 
//-------------
// C Headers --
//-------------
extern "C" {
#include <assert.h>
#include <math.h> 
#include <stdlib.h>
}
 
//---------------
// C++ Headers --
//---------------
#include <iostream>
//#include "ErrLogger/ErrLog.hh"
using std::endl;
 
double NumRecipes::gammln(double xx)
{
   double x,y,tmp,ser;
   static double cof[6]={76.18009172947146,-86.50532032941677,
       24.01409824083091,-1.231739572450155,
       0.1208650973866179e-2,-0.5395239384953e-5};
   int j;
 
   y=x=xx;
   tmp=x+5.5;
   tmp -= (x+0.5)*log(tmp);
   ser=1.000000000190015;
   for (j=0;j<=5;j++) ser += cof[j]/++y;
   return -tmp+log(2.5066282746310005*ser/x);
}
 
double 
NumRecipes::gammp(double a, double x)
{
   double gamser,gammcf,gln;
 
   if (x < 0.0 || a <= 0.0) std::cout<<"ErrMsg(error)" <<" Invalid arguments in routine gammp x=" << x << " a=" << a << endl;
   if (x < (a+1.0)) {
      gser(&gamser,a,x,&gln);
      return gamser;
   } else {
      gcf(&gammcf,a,x,&gln);
      return 1.0-gammcf;
   }
}
 
double 
NumRecipes::gammq(double a, double x)
{
   double gamser,gammcf,gln;
 
   if (x < 0.0 || a <= 0.0) recipesErr(" Invalid arguments in routine GAMMQ");
   if (x < (a+1.0)) {
      gser(&gamser,a,x,&gln);
      return 1.0-gamser;
   } else {
      gcf(&gammcf,a,x,&gln);
      return gammcf;
   }
}
 
void NumRecipes::gcf(double* gammcf, double a, double x, double* gln)
   {
     int n;
     double gold=0.0,g,fac=1.0,b1=1.0;
     double b0=0.0,anf,ana,an,a1,a0=1.0;
 
     *gln=gammln(a);
     a1=x;
     for (n=1;n<=NUMREC_ITMAX;n++) {
       an=(double) n;
       ana=an-a;
       a0=(a1+a0*ana)*fac;
       b0=(b1+b0*ana)*fac;
       anf=an*fac;
       a1=x*a0+anf*a1;
       b1=x*b0+anf*b1;
       if (a1) {
         fac=1.0/a1;
         g=b1*fac;
         if (fabs((g-gold)/g) < NUMREC_EPS) {
            *gammcf=exp(-x+a*log(x)-(*gln))*g;
            return;
         }
         gold=g;
       }
     }
     recipesErr(" a too large, NUMREC_ITMAX too small in routine GCF");
   }
 
void NumRecipes::gser(double* gamser, double a, double x, double* gln)
{
   int n;
   double sum,del,ap;
 
   *gln=gammln(a);
   if (x <= 0.0) {
      if (x < 0.0) recipesErr(" x less than 0 in routine GSER");
      *gamser=0.0;
      return;
   } else {
      ap=a;
      del=sum=1.0/a;
      for (n=1;n<=NUMREC_ITMAX;n++) {
         ap += 1.0;
         del *= x/ap;
         sum += del;
         if (fabs(del) < fabs(sum)*NUMREC_EPS) {
            *gamser=sum*exp(-x+a*log(x)-(*gln));
            return;
         }
      }
      recipesErr(" a too large, NUMREC_ITMAX too small in routine GSER");
      return;
   }
}
 
void NumRecipes::recipesErr(const char* c)
{
   std::cout<<"ErrMsg(fatal)" << " Numerical Recipes run-time error...\n" << c 
                 << "\n ...now exiting to system..." << std::endl;
   ::abort();
}