Ruby 1.9.3p327(2012-11-10revision37606)
|
00001 /* erf.c - public domain implementation of error function erf(3m) 00002 00003 reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten 00004 (New Algorithm handbook in C language) (Gijyutsu hyouron 00005 sha, Tokyo, 1991) p.227 [in Japanese] */ 00006 #include "ruby/missing.h" 00007 #include <stdio.h> 00008 #include <math.h> 00009 00010 #ifdef _WIN32 00011 # include <float.h> 00012 # if !defined __MINGW32__ || defined __NO_ISOCEXT 00013 # ifndef isnan 00014 # define isnan(x) _isnan(x) 00015 # endif 00016 # ifndef isinf 00017 # define isinf(x) (!_finite(x) && !_isnan(x)) 00018 # endif 00019 # ifndef finite 00020 # define finite(x) _finite(x) 00021 # endif 00022 # endif 00023 #endif 00024 00025 static double q_gamma(double, double, double); 00026 00027 /* Incomplete gamma function 00028 1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt */ 00029 static double p_gamma(double a, double x, double loggamma_a) 00030 { 00031 int k; 00032 double result, term, previous; 00033 00034 if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a); 00035 if (x == 0) return 0; 00036 result = term = exp(a * log(x) - x - loggamma_a) / a; 00037 for (k = 1; k < 1000; k++) { 00038 term *= x / (a + k); 00039 previous = result; result += term; 00040 if (result == previous) return result; 00041 } 00042 fprintf(stderr, "erf.c:%d:p_gamma() could not converge.", __LINE__); 00043 return result; 00044 } 00045 00046 /* Incomplete gamma function 00047 1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt */ 00048 static double q_gamma(double a, double x, double loggamma_a) 00049 { 00050 int k; 00051 double result, w, temp, previous; 00052 double la = 1, lb = 1 + x - a; /* Laguerre polynomial */ 00053 00054 if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a); 00055 w = exp(a * log(x) - x - loggamma_a); 00056 result = w / lb; 00057 for (k = 2; k < 1000; k++) { 00058 temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k; 00059 la = lb; lb = temp; 00060 w *= (k - 1 - a) / k; 00061 temp = w / (la * lb); 00062 previous = result; result += temp; 00063 if (result == previous) return result; 00064 } 00065 fprintf(stderr, "erf.c:%d:q_gamma() could not converge.", __LINE__); 00066 return result; 00067 } 00068 00069 #define LOG_PI_OVER_2 0.572364942924700087071713675675 /* log_e(PI)/2 */ 00070 00071 double erf(double x) 00072 { 00073 if (!finite(x)) { 00074 if (isnan(x)) return x; /* erf(NaN) = NaN */ 00075 return (x>0 ? 1.0 : -1.0); /* erf(+-inf) = +-1.0 */ 00076 } 00077 if (x >= 0) return p_gamma(0.5, x * x, LOG_PI_OVER_2); 00078 else return - p_gamma(0.5, x * x, LOG_PI_OVER_2); 00079 } 00080 00081 double erfc(double x) 00082 { 00083 if (!finite(x)) { 00084 if (isnan(x)) return x; /* erfc(NaN) = NaN */ 00085 return (x>0 ? 0.0 : 2.0); /* erfc(+-inf) = 0.0, 2.0 */ 00086 } 00087 if (x >= 0) return q_gamma(0.5, x * x, LOG_PI_OVER_2); 00088 else return 1 + p_gamma(0.5, x * x, LOG_PI_OVER_2); 00089 } 00090