libgfortran/intrinsics/erfc_scaled_inc.c

   1 /* Implementation of the ERFC_SCALED intrinsic, to be included by erfc_scaled.c
   2    Copyright (c) 2008 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 2 of the License, or (at your option) any later version.
  10
  11 In addition to the permissions in the GNU General Public License, the
  12 Free Software Foundation gives you unlimited permission to link the
  13 compiled version of this file into combinations with other programs,
  14 and to distribute those combinations without any restriction coming
  15 from the use of this file.  (The General Public License restrictions
  16 do apply in other respects; for example, they cover modification of
  17 the file, and distribution when not linked into a combine
  18 executable.)
  19
  20 Libgfortran is distributed in the hope that it will be useful,
  21 but WITHOUT ANY WARRANTY; without even the implied warranty of
  22 MERCHANTABILITY or FITNESS FOR a PARTICULAR PURPOSE.  See the
  23 GNU General Public License for more details.
  24
  25 You should have received a copy of the GNU General Public
  26 License along with libgfortran; see the file COPYING.  If not,
  27 write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  28 Boston, MA 02110-1301, USA.  */
  29
  30 /* This implementation of ERFC_SCALED is based on the netlib algorithm
  31    available at http://www.netlib.org/specfun/erf  */
  32
  33 #define TYPE KIND_SUFFIX(GFC_REAL_,KIND)
  34 #define CONCAT(x,y) x ## y
  35 #define KIND_SUFFIX(x,y) CONCAT(x,y)
  36
  37 #if (KIND == 4)
  38
  39 # define EXP(x) expf(x)
  40 # define TRUNC(x) truncf(x)
  41
  42 #elif (KIND == 8)
  43
  44 # define EXP(x) exp(x)
  45 # define TRUNC(x) trunc(x)
  46
  47 #else
  48
  49 # ifdef HAVE_EXPL
  50 #  define EXP(x) expl(x)
  51 # endif
  52 # ifdef HAVE_TRUNCL
  53 #  define TRUNC(x) truncl(x)
  54 # endif
  55
  56 #endif
  57
  58 #if defined(EXP) && defined(TRUNC)
  59
  60 extern TYPE KIND_SUFFIX(erfc_scaled_r,KIND) (TYPE);
  61 export_proto(KIND_SUFFIX(erfc_scaled_r,KIND));
  62
  63 TYPE
  64 KIND_SUFFIX(erfc_scaled_r,KIND) (TYPE x)
  65 {
  66   /* The main computation evaluates near-minimax approximations
  67      from "Rational Chebyshev approximations for the error function"
  68      by W. J. Cody, Math. Comp., 1969, PP. 631-638.  This
  69      transportable program uses rational functions that theoretically
  70      approximate  erf(x)  and  erfc(x)  to at least 18 significant
  71      decimal digits.  The accuracy achieved depends on the arithmetic
  72      system, the compiler, the intrinsic functions, and proper
  73      selection of the machine-dependent constants.  */
  74
  75   int i;
  76   TYPE del, res, xden, xnum, y, ysq;
  77
  78 #if (KIND == 4)
  79   static TYPE xneg = -9.382, xsmall = 5.96e-8,
  80               xbig = 9.194, xhuge = 2.90e+3, xmax = 4.79e+37;
  81 #else
  82   static TYPE xneg = -26.628, xsmall = 1.11e-16,
  83               xbig = 26.543, xhuge = 6.71e+7, xmax = 2.53e+307;
  84 #endif
  85
  86 #define SQRPI ((TYPE) 0.56418958354775628695L)
  87 #define THRESH ((TYPE) 0.46875L)
  88
  89   static TYPE a[5] = { 3.16112374387056560l, 113.864154151050156l,
  90     377.485237685302021l, 3209.37758913846947l, 0.185777706184603153l };
  91
  92   static TYPE b[4] = { 23.6012909523441209l, 244.024637934444173l,
  93     1282.61652607737228l, 2844.23683343917062l };
  94
  95   static TYPE c[9] = { 0.564188496988670089l, 8.88314979438837594l,
  96     66.1191906371416295l, 298.635138197400131l, 881.952221241769090l,
  97     1712.04761263407058l, 2051.07837782607147l, 1230.33935479799725l,
  98     2.15311535474403846e-8l };
  99
 100   static TYPE d[8] = { 15.7449261107098347l, 117.693950891312499l,
 101     537.181101862009858l, 1621.38957456669019l, 3290.79923573345963l,
 102     4362.61909014324716l, 3439.36767414372164l, 1230.33935480374942l };
 103
 104   static TYPE p[6] = { 0.305326634961232344l, 0.360344899949804439l,
 105     0.125781726111229246l, 0.0160837851487422766l,
 106     0.000658749161529837803l, 0.0163153871373020978l };
 107
 108   static TYPE q[5] = { 2.56852019228982242l, 1.87295284992346047l,
 109     0.527905102951428412l, 0.0605183413124413191l,
 110     0.00233520497626869185l };
 111
 112   y = (x > 0 ? x : -x);
 113   if (y <= THRESH)
 114     {
 115       ysq = 0;
 116       if (y > xsmall)
 117         ysq = y * y;
 118       xnum = a[4]*ysq;
 119       xden = ysq;
 120       for (i = 0; i <= 2; i++)
 121         {
 122           xnum = (xnum + a[i]) * ysq;
 123           xden = (xden + b[i]) * ysq;
 124         }
 125       res = x * (xnum + a[3]) / (xden + b[3]);
 126       res = 1 - res;
 127       res = EXP(ysq) * res;
 128       return res;
 129     }
 130   else if (y <= 4)
 131     {
 132       xnum = c[8]*y;
 133       xden = y;
 134       for (i = 0; i <= 6; i++)
 135         {
 136           xnum = (xnum + c[i]) * y;
 137           xden = (xden + d[i]) * y;
 138         }
 139       res = (xnum + c[7]) / (xden + d[7]);
 140     }
 141   else
 142     {
 143       res = 0;
 144       if (y >= xbig)
 145         {
 146           if (y >= xmax)
 147             goto finish;
 148           if (y >= xhuge)
 149             {
 150               res = SQRPI / y;
 151               goto finish;
 152             }
 153         }
 154       ysq = ((TYPE) 1) / (y * y);
 155       xnum = p[5]*ysq;
 156       xden = ysq;
 157       for (i = 0; i <= 3; i++)
 158         {
 159           xnum = (xnum + p[i]) * ysq;
 160           xden = (xden + q[i]) * ysq;
 161         }
 162       res = ysq *(xnum + p[4]) / (xden + q[4]);
 163       res = (SQRPI -  res) / y;
 164     }
 165
 166 finish:
 167   if (x < 0)
 168     {
 169       if (x < xneg)
 170         res = __builtin_inf ();
 171       else
 172         {
 173           ysq = TRUNC (x*((TYPE) 16))/((TYPE) 16);
 174           del = (x-ysq)*(x+ysq);
 175           y = EXP(ysq*ysq) * EXP(del);
 176           res = (y+y) - res;
 177         }
 178     }
 179   return res;
 180 }
 181
 182 #endif
 183
 184 #undef EXP
 185 #undef TRUNC
 186
 187 #undef CONCAT
 188 #undef TYPE
 189 #undef KIND_SUFFIX