1/* Implementation of the ERFC_SCALED intrinsic. 2 Copyright (C) 2008-2020 Free Software Foundation, Inc. 3 4This file is part of the GNU Fortran runtime library (libgfortran). 5 6Libgfortran is free software; you can redistribute it and/or 7modify it under the terms of the GNU General Public 8License as published by the Free Software Foundation; either 9version 3 of the License, or (at your option) any later version. 10 11Libgfortran is distributed in the hope that it will be useful, 12but WITHOUT ANY WARRANTY; without even the implied warranty of 13MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14GNU General Public License for more details. 15 16Under Section 7 of GPL version 3, you are granted additional 17permissions described in the GCC Runtime Library Exception, version 183.1, as published by the Free Software Foundation. 19 20You should have received a copy of the GNU General Public License and 21a copy of the GCC Runtime Library Exception along with this program; 22see 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 80extern GFC_REAL_16 erfc_scaled_r16 (GFC_REAL_16); 81export_proto(erfc_scaled_r16); 82 83GFC_REAL_16 84erfc_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