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