1 /* Implementation of the ERFC_SCALED intrinsic. 2 Copyright (C) 2008-2022 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 #define ERFC_SCALED(k) \ 79 GFC_REAL_ ## k \ 80 erfc_scaled_r ## k (GFC_REAL_ ## k x) \ 81 { \ 82 if (x < _THRESH) \ 83 { \ 84 return _INF; \ 85 } \ 86 if (x < 12) \ 87 { \ 88 /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2). \ 89 This is not perfect, but much better than netlib. */ \ 90 return _ERFC(x) * _EXP(x * x); \ 91 } \ 92 else \ 93 { \ 94 /* Calculate ERFC_SCALED(x) using a power series in 1/x: \ 95 ERFC_SCALED(x) = 1 / (x * sqrt(pi)) \ 96 * (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1)) \ 97 / (2 * x**2)**n) \ 98 */ \ 99 GFC_REAL_ ## k sum = 0, oldsum; \ 100 GFC_REAL_ ## k inv2x2 = 1 / (2 * x * x); \ 101 GFC_REAL_ ## k fac = 1; \ 102 int n = 1; \ 103 \ 104 while (n < 200) \ 105 { \ 106 fac *= - (2*n - 1) * inv2x2; \ 107 oldsum = sum; \ 108 sum += fac; \ 109 \ 110 if (sum == oldsum) \ 111 break; \ 112 \ 113 n++; \ 114 } \ 115 \ 116 return (1 + sum) / x * (_M_2_SQRTPI / 2); \ 117 } \ 118 } 119 120 #if defined(_ERFC) && defined(_EXP) 121 122 extern GFC_REAL_16 erfc_scaled_r16 (GFC_REAL_16); 123 export_proto(erfc_scaled_r16); 124 125 ERFC_SCALED(16) 126 127 #endif 128 129 #undef _THRESH 130 #undef _M_2_SQRTPI 131 #undef _INF 132 #undef _ERFC 133 #undef _EXP 134 135 #endif 136 137 #ifdef HAVE_GFC_REAL_17 138 139 /* For quadruple-precision, netlib's implementation is 140 not accurate enough. We provide another one. */ 141 142 # define _THRESH -106.566990228185312813205074546585730Q 143 # define _M_2_SQRTPI M_2_SQRTPIq 144 # define _INF __builtin_inff128() 145 # ifdef POWER_IEEE128 146 # define _ERFC(x) __erfcieee128(x) 147 # define _EXP(x) __expieee128(x) 148 # else 149 # define _ERFC(x) erfcq(x) 150 # define _EXP(x) expq(x) 151 # endif 152 153 extern GFC_REAL_17 erfc_scaled_r17 (GFC_REAL_17); 154 export_proto(erfc_scaled_r17); 155 156 ERFC_SCALED(17) 157 158 #undef _THRESH 159 #undef _M_2_SQRTPI 160 #undef _INF 161 #undef _ERFC 162 #undef _EXP 163 #undef ERFC_SCALED 164 165 #endif 166