summaryrefslogtreecommitdiff log msg author committer range
path: root/missing/erf.c
blob: d72c4eaf4e117cdfcbc57f2bc5036895a4c04dd2 (plain)
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 ``` ``````/* erf.c - public domain implementation of error function erf(3m) reference - Haruhiko Okumura: C-gengo niyoru saishin algorithm jiten (New Algorithm handbook in C language) (Gijyutsu hyouron sha, Tokyo, 1991) p.227 [in Japanese] */ #include "ruby/missing.h" #include #include #ifdef _WIN32 # include # if !defined __MINGW32__ || defined __NO_ISOCEXT # ifndef isnan # define isnan(x) _isnan(x) # endif # ifndef isinf # define isinf(x) (!_finite(x) && !_isnan(x)) # endif # ifndef finite # define finite(x) _finite(x) # endif # endif #endif static double q_gamma(double, double, double); /* Incomplete gamma function 1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt */ static double p_gamma(double a, double x, double loggamma_a) { int k; double result, term, previous; if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a); if (x == 0) return 0; result = term = exp(a * log(x) - x - loggamma_a) / a; for (k = 1; k < 1000; k++) { term *= x / (a + k); previous = result; result += term; if (result == previous) return result; } fprintf(stderr, "erf.c:%d:p_gamma() could not converge.", __LINE__); return result; } /* Incomplete gamma function 1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt */ static double q_gamma(double a, double x, double loggamma_a) { int k; double result, w, temp, previous; double la = 1, lb = 1 + x - a; /* Laguerre polynomial */ if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a); w = exp(a * log(x) - x - loggamma_a); result = w / lb; for (k = 2; k < 1000; k++) { temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k; la = lb; lb = temp; w *= (k - 1 - a) / k; temp = w / (la * lb); previous = result; result += temp; if (result == previous) return result; } fprintf(stderr, "erf.c:%d:q_gamma() could not converge.", __LINE__); return result; } #define LOG_PI_OVER_2 0.572364942924700087071713675675 /* log_e(PI)/2 */ double erf(double x) { if (!finite(x)) { if (isnan(x)) return x; /* erf(NaN) = NaN */ return (x>0 ? 1.0 : -1.0); /* erf(+-inf) = +-1.0 */ } if (x >= 0) return p_gamma(0.5, x * x, LOG_PI_OVER_2); else return - p_gamma(0.5, x * x, LOG_PI_OVER_2); } double erfc(double x) { if (!finite(x)) { if (isnan(x)) return x; /* erfc(NaN) = NaN */ return (x>0 ? 0.0 : 2.0); /* erfc(+-inf) = 0.0, 2.0 */ } if (x >= 0) return q_gamma(0.5, x * x, LOG_PI_OVER_2); else return 1 + p_gamma(0.5, x * x, LOG_PI_OVER_2); } ``````