1 /* Implementation of the ERFC_SCALED intrinsic, to be included by erfc_scaled.c
2 Copyright (C) 2008-2023 Free Software Foundation, Inc.
3
4 This file is part of the GNU Fortran runtime library (libgfortran).
5
6 Libgfortran is free software; you can redistribute it and/or
7 modify it under the terms of the GNU General Public
8 License as published by the Free Software Foundation; either
9 version 3 of the License, or (at your option) any later version.
10
11 Libgfortran is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR a PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 Under Section 7 of GPL version 3, you are granted additional
17 permissions described in the GCC Runtime Library Exception, version
18 3.1, as published by the Free Software Foundation.
19
20 You should have received a copy of the GNU General Public License and
21 a copy of the GCC Runtime Library Exception along with this program;
22 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
23 <http://www.gnu.org/licenses/>. */
24
25 /* This implementation of ERFC_SCALED is based on the netlib algorithm
26 available at http://www.netlib.org/specfun/erf */
27
28 #define TYPE KIND_SUFFIX(GFC_REAL_,KIND)
29 #define CONCAT(x,y) x ## y
30 #define KIND_SUFFIX(x,y) CONCAT(x,y)
31
32 #if (KIND == 4)
33
34 # define EXP(x) expf(x)
35 # define TRUNC(x) truncf(x)
36
37 #elif (KIND == 8)
38
39 # define EXP(x) exp(x)
40 # define TRUNC(x) trunc(x)
41
42 #elif (KIND == 10)
43
44 # ifdef HAVE_EXPL
45 # define EXP(x) expl(x)
46 # endif
47 # ifdef HAVE_TRUNCL
48 # define TRUNC(x) truncl(x)
49 # endif
50
51 #else
52
53 # error "What exactly is it that you want me to do?"
54
55 #endif
56
57 #if defined(EXP) && defined(TRUNC)
58
59 extern TYPE KIND_SUFFIX(erfc_scaled_r,KIND) (TYPE);
60 export_proto(KIND_SUFFIX(erfc_scaled_r,KIND));
61
62 TYPE
63 KIND_SUFFIX(erfc_scaled_r,KIND) (TYPE x)
64 {
65 /* The main computation evaluates near-minimax approximations
66 from "Rational Chebyshev approximations for the error function"
67 by W. J. Cody, Math. Comp., 1969, PP. 631-638. This
68 transportable program uses rational functions that theoretically
69 approximate erf(x) and erfc(x) to at least 18 significant
70 decimal digits. The accuracy achieved depends on the arithmetic
71 system, the compiler, the intrinsic functions, and proper
72 selection of the machine-dependent constants. */
73
74 int i;
75 TYPE del, res, xden, xnum, y, ysq;
76
77 #if (KIND == 4)
78 static TYPE xneg = -9.382, xsmall = 5.96e-8,
79 xbig = 9.194, xhuge = 2.90e+3, xmax = 4.79e+37;
80 #else
81 static TYPE xneg = -26.628, xsmall = 1.11e-16,
82 xbig = 26.543, xhuge = 6.71e+7, xmax = 2.53e+307;
83 #endif
84
85 #define SQRPI ((TYPE) 0.56418958354775628695L)
86 #define THRESH ((TYPE) 0.46875L)
87
88 static TYPE a[5] = { 3.16112374387056560l, 113.864154151050156l,
89 377.485237685302021l, 3209.37758913846947l, 0.185777706184603153l };
90
91 static TYPE b[4] = { 23.6012909523441209l, 244.024637934444173l,
92 1282.61652607737228l, 2844.23683343917062l };
93
94 static TYPE c[9] = { 0.564188496988670089l, 8.88314979438837594l,
95 66.1191906371416295l, 298.635138197400131l, 881.952221241769090l,
96 1712.04761263407058l, 2051.07837782607147l, 1230.33935479799725l,
97 2.15311535474403846e-8l };
98
99 static TYPE d[8] = { 15.7449261107098347l, 117.693950891312499l,
100 537.181101862009858l, 1621.38957456669019l, 3290.79923573345963l,
101 4362.61909014324716l, 3439.36767414372164l, 1230.33935480374942l };
102
103 static TYPE p[6] = { 0.305326634961232344l, 0.360344899949804439l,
104 0.125781726111229246l, 0.0160837851487422766l,
105 0.000658749161529837803l, 0.0163153871373020978l };
106
107 static TYPE q[5] = { 2.56852019228982242l, 1.87295284992346047l,
108 0.527905102951428412l, 0.0605183413124413191l,
109 0.00233520497626869185l };
110
111 y = (x > 0 ? x : -x);
112 if (y <= THRESH)
113 {
114 ysq = 0;
115 if (y > xsmall)
116 ysq = y * y;
117 xnum = a[4]*ysq;
118 xden = ysq;
119 for (i = 0; i <= 2; i++)
120 {
121 xnum = (xnum + a[i]) * ysq;
122 xden = (xden + b[i]) * ysq;
123 }
124 res = x * (xnum + a[3]) / (xden + b[3]);
125 res = 1 - res;
126 res = EXP(ysq) * res;
127 return res;
128 }
129 else if (y <= 4)
130 {
131 xnum = c[8]*y;
132 xden = y;
133 for (i = 0; i <= 6; i++)
134 {
135 xnum = (xnum + c[i]) * y;
136 xden = (xden + d[i]) * y;
137 }
138 res = (xnum + c[7]) / (xden + d[7]);
139 }
140 else
141 {
142 res = 0;
143 if (y >= xbig)
144 {
145 if (y >= xmax)
146 goto finish;
147 if (y >= xhuge)
148 {
149 res = SQRPI / y;
150 goto finish;
151 }
152 }
153 ysq = ((TYPE) 1) / (y * y);
154 xnum = p[5]*ysq;
155 xden = ysq;
156 for (i = 0; i <= 3; i++)
157 {
158 xnum = (xnum + p[i]) * ysq;
159 xden = (xden + q[i]) * ysq;
160 }
161 res = ysq *(xnum + p[4]) / (xden + q[4]);
162 res = (SQRPI - res) / y;
163 }
164
165 finish:
166 if (x < 0)
167 {
168 if (x < xneg)
169 res = __builtin_inf ();
170 else
171 {
172 ysq = TRUNC (x*((TYPE) 16))/((TYPE) 16);
173 del = (x-ysq)*(x+ysq);
174 y = EXP(ysq*ysq) * EXP(del);
175 res = (y+y) - res;
176 }
177 }
178 return res;
179 }
180
181 #endif
182
183 #undef EXP
184 #undef TRUNC
185
186 #undef CONCAT
187 #undef TYPE
188 #undef KIND_SUFFIX