libgfortran/intrinsics/erfc_scaled.c

   1 /* Implementation of the ERFC_SCALED intrinsic.
   2    Copyright (C) 2008-2023 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 # ifdef GFC_REAL_16_USE_IEC_60559
  56 #  define _THRESH -106.566990228185312813205074546585730F128
  57 #  define _M_2_SQRTPI M_2_SQRTPIf128
  58 #  define _INF __builtin_inff128()
  59 #  define _ERFC(x) erfcf128(x)
  60 #  define _EXP(x) expf128(x)
  61 # else
  62 #  define _THRESH -106.566990228185312813205074546585730Q
  63 #  define _M_2_SQRTPI M_2_SQRTPIq
  64 #  define _INF __builtin_infq()
  65 #  define _ERFC(x) erfcq(x)
  66 #  define _EXP(x) expq(x)
  67 # endif
  68
  69 #else
  70
  71 # define _THRESH -106.566990228185312813205074546585730L
  72 # ifndef M_2_SQRTPIl
  73 #  define M_2_SQRTPIl 1.128379167095512573896158903121545172L
  74 # endif
  75 # define _M_2_SQRTPI M_2_SQRTPIl
  76 # define _INF __builtin_infl()
  77 # ifdef HAVE_ERFCL
  78 #  define _ERFC(x) erfcl(x)
  79 # endif
  80 # ifdef HAVE_EXPL
  81 #  define _EXP(x) expl(x)
  82 # endif
  83
  84 #endif
  85
  86 #define ERFC_SCALED(k) \
  87 GFC_REAL_ ## k                                                              \
  88 erfc_scaled_r ## k (GFC_REAL_ ## k x)                                       \
  89 {                                                                           \
  90   if (x < _THRESH)                                                          \
  91     {                                                                       \
  92       return _INF;                                                          \
  93     }                                                                       \
  94   if (x < 12)                                                               \
  95     {                                                                       \
  96       /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2).          \
  97          This is not perfect, but much better than netlib.  */              \
  98       return _ERFC(x) * _EXP(x * x);                                        \
  99     }                                                                       \
 100   else                                                                      \
 101     {                                                                       \
 102       /* Calculate ERFC_SCALED(x) using a power series in 1/x:              \
 103          ERFC_SCALED(x) = 1 / (x * sqrt(pi))                                \
 104                          * (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1))  \
 105                                               / (2 * x**2)**n)              \
 106        */                                                                   \
 107       GFC_REAL_ ## k sum = 0, oldsum;                                       \
 108       GFC_REAL_ ## k inv2x2 = 1 / (2 * x * x);                              \
 109       GFC_REAL_ ## k fac = 1;                                               \
 110       int n = 1;                                                            \
 111                                                                             \
 112       while (n < 200)                                                       \
 113         {                                                                   \
 114           fac *= - (2*n - 1) * inv2x2;                                      \
 115           oldsum = sum;                                                     \
 116           sum += fac;                                                       \
 117                                                                             \
 118           if (sum == oldsum)                                                \
 119             break;                                                          \
 120                                                                             \
 121           n++;                                                              \
 122         }                                                                   \
 123                                                                             \
 124       return (1 + sum) / x * (_M_2_SQRTPI / 2);                             \
 125     }                                                                       \
 126 }
 127
 128 #if defined(_ERFC) && defined(_EXP)
 129
 130 extern GFC_REAL_16 erfc_scaled_r16 (GFC_REAL_16);
 131 export_proto(erfc_scaled_r16);
 132
 133 ERFC_SCALED(16)
 134
 135 #endif
 136
 137 #undef _THRESH
 138 #undef _M_2_SQRTPI
 139 #undef _INF
 140 #undef _ERFC
 141 #undef _EXP
 142
 143 #endif
 144
 145 #ifdef HAVE_GFC_REAL_17
 146
 147 /* For quadruple-precision, netlib's implementation is
 148    not accurate enough.  We provide another one.  */
 149
 150 # define _THRESH GFC_REAL_17_LITERAL(-106.566990228185312813205074546585730)
 151 # define _M_2_SQRTPI GFC_REAL_17_LITERAL(M_2_SQRTPI)
 152 # define _INF __builtin_inff128()
 153 # ifdef POWER_IEEE128
 154 #  define _ERFC(x) __erfcieee128(x)
 155 #  define _EXP(x) __expieee128(x)
 156 # elif defined(GFC_REAL_17_USE_IEC_60559)
 157 #  define _ERFC(x) erfcf128(x)
 158 #  define _EXP(x) expf128(x)
 159 # else
 160 #  define _ERFC(x) erfcq(x)
 161 #  define _EXP(x) expq(x)
 162 # endif
 163
 164 extern GFC_REAL_17 erfc_scaled_r17 (GFC_REAL_17);
 165 export_proto(erfc_scaled_r17);
 166
 167 ERFC_SCALED(17)
 168
 169 #undef _THRESH
 170 #undef _M_2_SQRTPI
 171 #undef _INF
 172 #undef _ERFC
 173 #undef _EXP
 174 #undef ERFC_SCALED
 175
 176 #endif