]>
Commit | Line | Data |
---|---|---|
f489fba1 | 1 | /* Implementation of the ERFC_SCALED intrinsic, to be included by erfc_scaled.c |
a5544970 | 2 | Copyright (C) 2008-2019 Free Software Foundation, Inc. |
f489fba1 FXC |
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 | |
748086b7 | 9 | version 3 of the License, or (at your option) any later version. |
f489fba1 FXC |
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 | ||
748086b7 JJ |
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/>. */ | |
f489fba1 FXC |
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) | |
cb31c4bc | 33 | |
f489fba1 FXC |
34 | # define EXP(x) expf(x) |
35 | # define TRUNC(x) truncf(x) | |
cb31c4bc | 36 | |
f489fba1 | 37 | #elif (KIND == 8) |
cb31c4bc | 38 | |
f489fba1 FXC |
39 | # define EXP(x) exp(x) |
40 | # define TRUNC(x) trunc(x) | |
cb31c4bc | 41 | |
3d41d9d9 | 42 | #elif (KIND == 10) |
cb31c4bc FXC |
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 | ||
1ec601bf FXC |
51 | #else |
52 | ||
53 | # error "What exactly is it that you want me to do?" | |
54 | ||
f489fba1 FXC |
55 | #endif |
56 | ||
cb31c4bc FXC |
57 | #if defined(EXP) && defined(TRUNC) |
58 | ||
f489fba1 FXC |
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 | ||
cb31c4bc FXC |
181 | #endif |
182 | ||
f489fba1 FXC |
183 | #undef EXP |
184 | #undef TRUNC | |
185 | ||
186 | #undef CONCAT | |
187 | #undef TYPE | |
188 | #undef KIND_SUFFIX |