00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "ProbTools/NumRecipes.h"
00023
00024
00025
00026
00027 extern "C" {
00028 #include <assert.h>
00029 #include <math.h>
00030 #include <stdlib.h>
00031 }
00032
00033
00034
00035
00036 #include <iostream>
00037
00038 using std::endl;
00039
00040 double NumRecipes::gammln(double xx)
00041 {
00042 double x,y,tmp,ser;
00043 static double cof[6]={76.18009172947146,-86.50532032941677,
00044 24.01409824083091,-1.231739572450155,
00045 0.1208650973866179e-2,-0.5395239384953e-5};
00046 int j;
00047
00048 y=x=xx;
00049 tmp=x+5.5;
00050 tmp -= (x+0.5)*log(tmp);
00051 ser=1.000000000190015;
00052 for (j=0;j<=5;j++) ser += cof[j]/++y;
00053 return -tmp+log(2.5066282746310005*ser/x);
00054 }
00055
00056 double
00057 NumRecipes::gammp(double a, double x)
00058 {
00059 double gamser,gammcf,gln;
00060
00061 if (x < 0.0 || a <= 0.0) std::cout<<"ErrMsg(error)" <<" Invalid arguments in routine gammp x=" << x << " a=" << a << endl;
00062 if (x < (a+1.0)) {
00063 gser(&gamser,a,x,&gln);
00064 return gamser;
00065 } else {
00066 gcf(&gammcf,a,x,&gln);
00067 return 1.0-gammcf;
00068 }
00069 }
00070
00071 double
00072 NumRecipes::gammq(double a, double x)
00073 {
00074 double gamser,gammcf,gln;
00075
00076 if (x < 0.0 || a <= 0.0) recipesErr(" Invalid arguments in routine GAMMQ");
00077 if (x < (a+1.0)) {
00078 gser(&gamser,a,x,&gln);
00079 return 1.0-gamser;
00080 } else {
00081 gcf(&gammcf,a,x,&gln);
00082 return gammcf;
00083 }
00084 }
00085
00086 void NumRecipes::gcf(double* gammcf, double a, double x, double* gln)
00087 {
00088 int n;
00089 double gold=0.0,g,fac=1.0,b1=1.0;
00090 double b0=0.0,anf,ana,an,a1,a0=1.0;
00091
00092 *gln=gammln(a);
00093 a1=x;
00094 for (n=1;n<=NUMREC_ITMAX;n++) {
00095 an=(double) n;
00096 ana=an-a;
00097 a0=(a1+a0*ana)*fac;
00098 b0=(b1+b0*ana)*fac;
00099 anf=an*fac;
00100 a1=x*a0+anf*a1;
00101 b1=x*b0+anf*b1;
00102 if (a1) {
00103 fac=1.0/a1;
00104 g=b1*fac;
00105 if (fabs((g-gold)/g) < NUMREC_EPS) {
00106 *gammcf=exp(-x+a*log(x)-(*gln))*g;
00107 return;
00108 }
00109 gold=g;
00110 }
00111 }
00112 recipesErr(" a too large, NUMREC_ITMAX too small in routine GCF");
00113 }
00114
00115 void NumRecipes::gser(double* gamser, double a, double x, double* gln)
00116 {
00117 int n;
00118 double sum,del,ap;
00119
00120 *gln=gammln(a);
00121 if (x <= 0.0) {
00122 if (x < 0.0) recipesErr(" x less than 0 in routine GSER");
00123 *gamser=0.0;
00124 return;
00125 } else {
00126 ap=a;
00127 del=sum=1.0/a;
00128 for (n=1;n<=NUMREC_ITMAX;n++) {
00129 ap += 1.0;
00130 del *= x/ap;
00131 sum += del;
00132 if (fabs(del) < fabs(sum)*NUMREC_EPS) {
00133 *gamser=sum*exp(-x+a*log(x)-(*gln));
00134 return;
00135 }
00136 }
00137 recipesErr(" a too large, NUMREC_ITMAX too small in routine GSER");
00138 return;
00139 }
00140 }
00141
00142 void NumRecipes::recipesErr(const char* c)
00143 {
00144 std::cout<<"ErrMsg(fatal)" << " Numerical Recipes run-time error...\n" << c
00145 << "\n ...now exiting to system..." << std::endl;
00146 ::abort();
00147 }
00148