missing/erf.c

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