1 /* Implementation of the ERFC_SCALED intrinsic.
2 Copyright (C) 2008-2019 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 # define _THRESH -106.566990228185312813205074546585730Q
56 # define _M_2_SQRTPI M_2_SQRTPIq
57 # define _INF __builtin_infq()
58 # define _ERFC(x) erfcq(x)
59 # define _EXP(x) expq(x)
60
61 #else
62
63 # define _THRESH -106.566990228185312813205074546585730L
64 # ifndef M_2_SQRTPIl
65 # define M_2_SQRTPIl 1.128379167095512573896158903121545172L
66 # endif
67 # define _M_2_SQRTPI M_2_SQRTPIl
68 # define _INF __builtin_infl()
69 # ifdef HAVE_ERFCL
70 # define _ERFC(x) erfcl(x)
71 # endif
72 # ifdef HAVE_EXPL
73 # define _EXP(x) expl(x)
74 # endif
75
76 #endif
77
78 #if defined(_ERFC) && defined(_EXP)
79
80 extern GFC_REAL_16 erfc_scaled_r16 (GFC_REAL_16);
81 export_proto(erfc_scaled_r16);
82
83 GFC_REAL_16
erfc_scaled_r16(GFC_REAL_16 x)84 erfc_scaled_r16 (GFC_REAL_16 x)
85 {
86 if (x < _THRESH)
87 {
88 return _INF;
89 }
90 if (x < 12)
91 {
92 /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2).
93 This is not perfect, but much better than netlib. */
94 return _ERFC(x) * _EXP(x * x);
95 }
96 else
97 {
98 /* Calculate ERFC_SCALED(x) using a power series in 1/x:
99 ERFC_SCALED(x) = 1 / (x * sqrt(pi))
100 * (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1))
101 / (2 * x**2)**n)
102 */
103 GFC_REAL_16 sum = 0, oldsum;
104 GFC_REAL_16 inv2x2 = 1 / (2 * x * x);
105 GFC_REAL_16 fac = 1;
106 int n = 1;
107
108 while (n < 200)
109 {
110 fac *= - (2*n - 1) * inv2x2;
111 oldsum = sum;
112 sum += fac;
113
114 if (sum == oldsum)
115 break;
116
117 n++;
118 }
119
120 return (1 + sum) / x * (_M_2_SQRTPI / 2);
121 }
122 }
123
124 #endif
125
126 #endif
127