1 /* Implementation of the ERFC_SCALED intrinsic.
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 #include "libgfortran.h"
26
27 /* This implementation of ERFC_SCALED is based on the netlib algorithm
28 available at http://www.netlib.org/specfun/erf */
29
30 #ifdef HAVE_GFC_REAL_4
31 #undef KIND
32 #define KIND 4
33 #include "erfc_scaled_inc.c"
34 #endif
35
36 #ifdef HAVE_GFC_REAL_8
37 #undef KIND
38 #define KIND 8
39 #include "erfc_scaled_inc.c"
40 #endif
41
42 #ifdef HAVE_GFC_REAL_10
43 #undef KIND
44 #define KIND 10
45 #include "erfc_scaled_inc.c"
46 #endif
47
48 #ifdef HAVE_GFC_REAL_16
49
50 /* For quadruple-precision, netlib's implementation is
51 not accurate enough. We provide another one. */
52
53 #ifdef GFC_REAL_16_IS_FLOAT128
54
55 # ifdef GFC_REAL_16_USE_IEC_60559
56 # define _THRESH -106.566990228185312813205074546585730F128
57 # define _M_2_SQRTPI M_2_SQRTPIf128
58 # define _INF __builtin_inff128()
59 # define _ERFC(x) erfcf128(x)
60 # define _EXP(x) expf128(x)
61 # else
62 # define _THRESH -106.566990228185312813205074546585730Q
63 # define _M_2_SQRTPI M_2_SQRTPIq
64 # define _INF __builtin_infq()
65 # define _ERFC(x) erfcq(x)
66 # define _EXP(x) expq(x)
67 # endif
68
69 #else
70
71 # define _THRESH -106.566990228185312813205074546585730L
72 # ifndef M_2_SQRTPIl
73 # define M_2_SQRTPIl 1.128379167095512573896158903121545172L
74 # endif
75 # define _M_2_SQRTPI M_2_SQRTPIl
76 # define _INF __builtin_infl()
77 # ifdef HAVE_ERFCL
78 # define _ERFC(x) erfcl(x)
79 # endif
80 # ifdef HAVE_EXPL
81 # define _EXP(x) expl(x)
82 # endif
83
84 #endif
85
86 #define ERFC_SCALED(k) \
87 GFC_REAL_ ## k \
88 erfc_scaled_r ## k (GFC_REAL_ ## k x) \
89 { \
90 if (x < _THRESH) \
91 { \
92 return _INF; \
93 } \
94 if (x < 12) \
95 { \
96 /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2). \
97 This is not perfect, but much better than netlib. */ \
98 return _ERFC(x) * _EXP(x * x); \
99 } \
100 else \
101 { \
102 /* Calculate ERFC_SCALED(x) using a power series in 1/x: \
103 ERFC_SCALED(x) = 1 / (x * sqrt(pi)) \
104 * (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1)) \
105 / (2 * x**2)**n) \
106 */ \
107 GFC_REAL_ ## k sum = 0, oldsum; \
108 GFC_REAL_ ## k inv2x2 = 1 / (2 * x * x); \
109 GFC_REAL_ ## k fac = 1; \
110 int n = 1; \
111 \
112 while (n < 200) \
113 { \
114 fac *= - (2*n - 1) * inv2x2; \
115 oldsum = sum; \
116 sum += fac; \
117 \
118 if (sum == oldsum) \
119 break; \
120 \
121 n++; \
122 } \
123 \
124 return (1 + sum) / x * (_M_2_SQRTPI / 2); \
125 } \
126 }
127
128 #if defined(_ERFC) && defined(_EXP)
129
130 extern GFC_REAL_16 erfc_scaled_r16 (GFC_REAL_16);
131 export_proto(erfc_scaled_r16);
132
133 ERFC_SCALED(16)
134
135 #endif
136
137 #undef _THRESH
138 #undef _M_2_SQRTPI
139 #undef _INF
140 #undef _ERFC
141 #undef _EXP
142
143 #endif
144
145 #ifdef HAVE_GFC_REAL_17
146
147 /* For quadruple-precision, netlib's implementation is
148 not accurate enough. We provide another one. */
149
150 # define _THRESH GFC_REAL_17_LITERAL(-106.566990228185312813205074546585730)
151 # define _M_2_SQRTPI GFC_REAL_17_LITERAL(M_2_SQRTPI)
152 # define _INF __builtin_inff128()
153 # ifdef POWER_IEEE128
154 # define _ERFC(x) __erfcieee128(x)
155 # define _EXP(x) __expieee128(x)
156 # elif defined(GFC_REAL_17_USE_IEC_60559)
157 # define _ERFC(x) erfcf128(x)
158 # define _EXP(x) expf128(x)
159 # else
160 # define _ERFC(x) erfcq(x)
161 # define _EXP(x) expq(x)
162 # endif
163
164 extern GFC_REAL_17 erfc_scaled_r17 (GFC_REAL_17);
165 export_proto(erfc_scaled_r17);
166
167 ERFC_SCALED(17)
168
169 #undef _THRESH
170 #undef _M_2_SQRTPI
171 #undef _INF
172 #undef _ERFC
173 #undef _EXP
174 #undef ERFC_SCALED
175
176 #endif