libgfortran/intrinsics/erfc_scaled_inc.c

   1 /* Implementation of the ERFC_SCALED intrinsic, to be included by erfc_scaled.c
   2    Copyright (C) 2008-2014 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 /* This implementation of ERFC_SCALED is based on the netlib algorithm
  26    available at http://www.netlib.org/specfun/erf  */
  27
  28 #define TYPE KIND_SUFFIX(GFC_REAL_,KIND)
  29 #define CONCAT(x,y) x ## y
  30 #define KIND_SUFFIX(x,y) CONCAT(x,y)
  31
  32 #if (KIND == 4)
  33
  34 # define EXP(x) expf(x)
  35 # define TRUNC(x) truncf(x)
  36
  37 #elif (KIND == 8)
  38
  39 # define EXP(x) exp(x)
  40 # define TRUNC(x) trunc(x)
  41
  42 #elif (KIND == 10)
  43
  44 # ifdef HAVE_EXPL
  45 #  define EXP(x) expl(x)
  46 # endif
  47 # ifdef HAVE_TRUNCL
  48 #  define TRUNC(x) truncl(x)
  49 # endif
  50
  51 #else
  52
  53 # error "What exactly is it that you want me to do?"
  54
  55 #endif
  56
  57 #if defined(EXP) && defined(TRUNC)
  58
  59 extern TYPE KIND_SUFFIX(erfc_scaled_r,KIND) (TYPE);
  60 export_proto(KIND_SUFFIX(erfc_scaled_r,KIND));
  61
  62 TYPE
  63 KIND_SUFFIX(erfc_scaled_r,KIND) (TYPE x)
  64 {
  65   /* The main computation evaluates near-minimax approximations
  66      from "Rational Chebyshev approximations for the error function"
  67      by W. J. Cody, Math. Comp., 1969, PP. 631-638.  This
  68      transportable program uses rational functions that theoretically
  69      approximate  erf(x)  and  erfc(x)  to at least 18 significant
  70      decimal digits.  The accuracy achieved depends on the arithmetic
  71      system, the compiler, the intrinsic functions, and proper
  72      selection of the machine-dependent constants.  */
  73
  74   int i;
  75   TYPE del, res, xden, xnum, y, ysq;
  76
  77 #if (KIND == 4)
  78   static TYPE xneg = -9.382, xsmall = 5.96e-8,
  79               xbig = 9.194, xhuge = 2.90e+3, xmax = 4.79e+37;
  80 #else
  81   static TYPE xneg = -26.628, xsmall = 1.11e-16,
  82               xbig = 26.543, xhuge = 6.71e+7, xmax = 2.53e+307;
  83 #endif
  84
  85 #define SQRPI ((TYPE) 0.56418958354775628695L)
  86 #define THRESH ((TYPE) 0.46875L)
  87
  88   static TYPE a[5] = { 3.16112374387056560l, 113.864154151050156l,
  89     377.485237685302021l, 3209.37758913846947l, 0.185777706184603153l };
  90
  91   static TYPE b[4] = { 23.6012909523441209l, 244.024637934444173l,
  92     1282.61652607737228l, 2844.23683343917062l };
  93
  94   static TYPE c[9] = { 0.564188496988670089l, 8.88314979438837594l,
  95     66.1191906371416295l, 298.635138197400131l, 881.952221241769090l,
  96     1712.04761263407058l, 2051.07837782607147l, 1230.33935479799725l,
  97     2.15311535474403846e-8l };
  98
  99   static TYPE d[8] = { 15.7449261107098347l, 117.693950891312499l,
 100     537.181101862009858l, 1621.38957456669019l, 3290.79923573345963l,
 101     4362.61909014324716l, 3439.36767414372164l, 1230.33935480374942l };
 102
 103   static TYPE p[6] = { 0.305326634961232344l, 0.360344899949804439l,
 104     0.125781726111229246l, 0.0160837851487422766l,
 105     0.000658749161529837803l, 0.0163153871373020978l };
 106
 107   static TYPE q[5] = { 2.56852019228982242l, 1.87295284992346047l,
 108     0.527905102951428412l, 0.0605183413124413191l,
 109     0.00233520497626869185l };
 110
 111   y = (x > 0 ? x : -x);
 112   if (y <= THRESH)
 113     {
 114       ysq = 0;
 115       if (y > xsmall)
 116         ysq = y * y;
 117       xnum = a[4]*ysq;
 118       xden = ysq;
 119       for (i = 0; i <= 2; i++)
 120         {
 121           xnum = (xnum + a[i]) * ysq;
 122           xden = (xden + b[i]) * ysq;
 123         }
 124       res = x * (xnum + a[3]) / (xden + b[3]);
 125       res = 1 - res;
 126       res = EXP(ysq) * res;
 127       return res;
 128     }
 129   else if (y <= 4)
 130     {
 131       xnum = c[8]*y;
 132       xden = y;
 133       for (i = 0; i <= 6; i++)
 134         {
 135           xnum = (xnum + c[i]) * y;
 136           xden = (xden + d[i]) * y;
 137         }
 138       res = (xnum + c[7]) / (xden + d[7]);
 139     }
 140   else
 141     {
 142       res = 0;
 143       if (y >= xbig)
 144         {
 145           if (y >= xmax)
 146             goto finish;
 147           if (y >= xhuge)
 148             {
 149               res = SQRPI / y;
 150               goto finish;
 151             }
 152         }
 153       ysq = ((TYPE) 1) / (y * y);
 154       xnum = p[5]*ysq;
 155       xden = ysq;
 156       for (i = 0; i <= 3; i++)
 157         {
 158           xnum = (xnum + p[i]) * ysq;
 159           xden = (xden + q[i]) * ysq;
 160         }
 161       res = ysq *(xnum + p[4]) / (xden + q[4]);
 162       res = (SQRPI -  res) / y;
 163     }
 164
 165 finish:
 166   if (x < 0)
 167     {
 168       if (x < xneg)
 169         res = __builtin_inf ();
 170       else
 171         {
 172           ysq = TRUNC (x*((TYPE) 16))/((TYPE) 16);
 173           del = (x-ysq)*(x+ysq);
 174           y = EXP(ysq*ysq) * EXP(del);
 175           res = (y+y) - res;
 176         }
 177     }
 178   return res;
 179 }
 180
 181 #endif
 182
 183 #undef EXP
 184 #undef TRUNC
 185
 186 #undef CONCAT
 187 #undef TYPE
 188 #undef KIND_SUFFIX