libgfortran/intrinsics/erfc_scaled.c

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