Go to the documentation of this file.00001
00002
00003
00004
00005
00006 #include <stdio.h>
00007 #include <math.h>
00008
00009 #ifdef _WIN32
00010 # include <float.h>
00011 # if !defined __MINGW32__ || defined __NO_ISOCEXT
00012 # ifndef isnan
00013 # define isnan(x) _isnan(x)
00014 # endif
00015 # ifndef isinf
00016 # define isinf(x) (!_finite(x) && !_isnan(x))
00017 # endif
00018 # ifndef finite
00019 # define finite(x) _finite(x)
00020 # endif
00021 # endif
00022 #endif
00023
00024 static double q_gamma(double, double, double);
00025
00026
00027
00028 static double p_gamma(double a, double x, double loggamma_a)
00029 {
00030 int k;
00031 double result, term, previous;
00032
00033 if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);
00034 if (x == 0) return 0;
00035 result = term = exp(a * log(x) - x - loggamma_a) / a;
00036 for (k = 1; k < 1000; k++) {
00037 term *= x / (a + k);
00038 previous = result; result += term;
00039 if (result == previous) return result;
00040 }
00041 fprintf(stderr, "erf.c:%d:p_gamma() could not converge.", __LINE__);
00042 return result;
00043 }
00044
00045
00046
00047 static double q_gamma(double a, double x, double loggamma_a)
00048 {
00049 int k;
00050 double result, w, temp, previous;
00051 double la = 1, lb = 1 + x - a;
00052
00053 if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);
00054 w = exp(a * log(x) - x - loggamma_a);
00055 result = w / lb;
00056 for (k = 2; k < 1000; k++) {
00057 temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;
00058 la = lb; lb = temp;
00059 w *= (k - 1 - a) / k;
00060 temp = w / (la * lb);
00061 previous = result; result += temp;
00062 if (result == previous) return result;
00063 }
00064 fprintf(stderr, "erf.c:%d:q_gamma() could not converge.", __LINE__);
00065 return result;
00066 }
00067
00068 #define LOG_PI_OVER_2 0.572364942924700087071713675675
00069
00070 double erf(double x)
00071 {
00072 if (!finite(x)) {
00073 if (isnan(x)) return x;
00074 return (x>0 ? 1.0 : -1.0);
00075 }
00076 if (x >= 0) return p_gamma(0.5, x * x, LOG_PI_OVER_2);
00077 else return - p_gamma(0.5, x * x, LOG_PI_OVER_2);
00078 }
00079
00080 double erfc(double x)
00081 {
00082 if (!finite(x)) {
00083 if (isnan(x)) return x;
00084 return (x>0 ? 0.0 : 2.0);
00085 }
00086 if (x >= 0) return q_gamma(0.5, x * x, LOG_PI_OVER_2);
00087 else return 1 + p_gamma(0.5, x * x, LOG_PI_OVER_2);
00088 }
00089