Ruby  3.1.4p223 (2023-03-30 revision HEAD)
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 "ruby/missing.h"
7 #include <stdio.h>
8 #include <math.h>
9 
10 static double q_gamma(double, double, double);
11 
12 /* Incomplete gamma function
13  1 / Gamma(a) * Int_0^x exp(-t) t^(a-1) dt */
14 static double p_gamma(double a, double x, double loggamma_a)
15 {
16  int k;
17  double result, term, previous;
18 
19  if (x >= 1 + a) return 1 - q_gamma(a, x, loggamma_a);
20  if (x == 0) return 0;
21  result = term = exp(a * log(x) - x - loggamma_a) / a;
22  for (k = 1; k < 1000; k++) {
23  term *= x / (a + k);
24  previous = result; result += term;
25  if (result == previous) return result;
26  }
27  fprintf(stderr, "erf.c:%d:p_gamma() could not converge.", __LINE__);
28  return result;
29 }
30 
31 /* Incomplete gamma function
32  1 / Gamma(a) * Int_x^inf exp(-t) t^(a-1) dt */
33 static double q_gamma(double a, double x, double loggamma_a)
34 {
35  int k;
36  double result, w, temp, previous;
37  double la = 1, lb = 1 + x - a; /* Laguerre polynomial */
38 
39  if (x < 1 + a) return 1 - p_gamma(a, x, loggamma_a);
40  w = exp(a * log(x) - x - loggamma_a);
41  result = w / lb;
42  for (k = 2; k < 1000; k++) {
43  temp = ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;
44  la = lb; lb = temp;
45  w *= (k - 1 - a) / k;
46  temp = w / (la * lb);
47  previous = result; result += temp;
48  if (result == previous) return result;
49  }
50  fprintf(stderr, "erf.c:%d:q_gamma() could not converge.", __LINE__);
51  return result;
52 }
53 
54 #define LOG_PI_OVER_2 0.572364942924700087071713675675 /* log_e(PI)/2 */
55 
56 double erf(double x)
57 {
58  if (!finite(x)) {
59  if (isnan(x)) return x; /* erf(NaN) = NaN */
60  return (x>0 ? 1.0 : -1.0); /* erf(+-inf) = +-1.0 */
61  }
62  if (x >= 0) return p_gamma(0.5, x * x, LOG_PI_OVER_2);
63  else return - p_gamma(0.5, x * x, LOG_PI_OVER_2);
64 }
65 
66 double erfc(double x)
67 {
68  if (!finite(x)) {
69  if (isnan(x)) return x; /* erfc(NaN) = NaN */
70  return (x>0 ? 0.0 : 2.0); /* erfc(+-inf) = 0.0, 2.0 */
71  }
72  if (x >= 0) return q_gamma(0.5, x * x, LOG_PI_OVER_2);
73  else return 1 + p_gamma(0.5, x * x, LOG_PI_OVER_2);
74 }