00001
00008 #include <string>
00009 #include "erf.h"
00010 #include <iostream>
00011
00012 using namespace std;
00013
00014 bool str2bool(string str);
00015 string bool2str(bool b);
00016
00017 void nrerror(const char * msg) {
00018 cout << msg;
00019 return;
00020 }
00021
00023 double gammln(double xx)
00024 {
00025
00026 double x,y,tmp,ser;
00027 static double cof[6]={76.18009172947146,-86.50532032941677, 24.01409824083091,-1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};
00028 int j;
00029 y=x=xx;
00030 tmp=x+5.5;
00031 tmp -= (x+0.5)*log(tmp);
00032 ser=1.000000000190015;
00033 for (j=0;j<=5;j++) ser += cof[j]/++y;
00034 return -tmp+log(2.5066282746310005*ser/x);
00035 }
00036
00037
00039 double gammp(double a, double x)
00040 {
00041 void gcf(double *gammcf, double a, double x, double *gln);
00042 void gser(double *gamser, double a, double x, double *gln);
00043
00044 double gamser,gammcf,gln;
00045 if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammp");
00046 if (x < (a+1.0)) {
00047 gser(&gamser,a,x,&gln);
00048 return gamser;
00049 } else {
00050 gcf(&gammcf,a,x,&gln);
00051 return 1.0-gammcf;
00052 }
00053 }
00054
00056 double gammq(double a, double x)
00057 {
00058 void gcf(double *gammcf, double a, double x, double *gln);
00059 void gser(double *gamser, double a, double x, double *gln);
00060
00061 double gamser,gammcf,gln;
00062 if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammq");
00063 if (x < (a+1.0)) {
00064 gser(&gamser,a,x,&gln);
00065 return 1.0-gamser;
00066 } else {
00067 gcf(&gammcf,a,x,&gln);
00068 return gammcf;
00069 }
00070 }
00071
00074 void gser(double *gamser, double a, double x, double *gln)
00075 {
00076 double gammln(double xx);
00077
00078 int n;
00079 double sum,del,ap;
00080 *gln=gammln(a);
00081 if (x <= 0.0) {
00082 if (x < 0.0) nrerror("x less than 0 in routine gser");
00083 *gamser=0.0;
00084 return;
00085 } else {
00086 ap=a;
00087 del=sum=1.0/a;
00088 for (n=1;n<=ITMAX;n++) {
00089 ++ap;
00090 del *= x/ap;
00091 sum += del;
00092 if (fabs(del) < fabs(sum)*EPS) {
00093 *gamser=sum*exp(-x+a*log(x)-(*gln));
00094 return;
00095 }
00096 }
00097 nrerror("a too large, ITMAX too small in routine gser");
00098 return;
00099 }
00100 }
00101
00103 void gcf(double *gammcf, double a, double x, double *gln)
00104 {
00105 double gammln(double xx);
00106 double an;
00107 double b;
00108 double c1;
00109 double d;
00110 double del;
00111 double h;
00112
00113 int i;
00114 *gln=gammln(a);
00115 b=x+1.0-a;
00116 c1=1.0/FPMIN;
00117 d=1.0/b;
00118 h=d;
00119 for (i=1;i<=ITMAX;i++) {
00120 an = -i*(i-a);
00121 b += 2.0;
00122 d=an*d+b;
00123 if (fabs(d) < FPMIN) d=FPMIN;
00124 c1=b+an/c1;
00125 if (fabs(c1) < FPMIN) c1=FPMIN;
00126 d=1.0/d;
00127 del=d*c1;
00128 h *= del;
00129 if (fabs(del-1.0) < EPS) break;
00130 }
00131 if (i > ITMAX) nrerror("a too large, ITMAX too small in gcf");
00132 *gammcf=exp(-x+a*log(x)-(*gln))*h;
00133 }
00134
00135
00137 double erf(double x)
00138 {
00139 double gammp(double a, double x);
00140 return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
00141 }