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