|
111
|
1 /* Implementation of the ERFC_SCALED intrinsic.
|
|
131
|
2 Copyright (C) 2008-2018 Free Software Foundation, Inc.
|
|
111
|
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
|
|
|
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
|
