Ruby  2.0.0p247(2013-06-27revision41674)
missing/erf.c
Go to the documentation of this file.
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