libgfortran/intrinsics/erfc_scaled_inc.c

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