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