1 /* Implementation of the ERFC_SCALED intrinsic. 2 Copyright (C) 2008-2020 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 127