Mercurial > hg > CbC > CbC_gcc
comparison libgfortran/intrinsics/erfc_scaled.c @ 111:04ced10e8804
gcc 7
author | kono |
---|---|
date | Fri, 27 Oct 2017 22:46:09 +0900 |
parents | |
children | 84e7813d76e9 |
comparison
equal
deleted
inserted
replaced
68:561a7518be6b | 111:04ced10e8804 |
---|---|
1 /* Implementation of the ERFC_SCALED intrinsic. | |
2 Copyright (C) 2008-2017 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 | |
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 |