diff options
Diffstat (limited to 'libgfortran/intrinsics/erfc_scaled_inc.c')
| -rw-r--r-- | libgfortran/intrinsics/erfc_scaled_inc.c | 193 |
1 files changed, 193 insertions, 0 deletions
diff --git a/libgfortran/intrinsics/erfc_scaled_inc.c b/libgfortran/intrinsics/erfc_scaled_inc.c new file mode 100644 index 000000000..7886136c5 --- /dev/null +++ b/libgfortran/intrinsics/erfc_scaled_inc.c @@ -0,0 +1,193 @@ +/* Implementation of the ERFC_SCALED intrinsic, to be included by erfc_scaled.c + Copyright (c) 2008, 2010 Free Software Foundation, Inc. + +This file is part of the GNU Fortran runtime library (libgfortran). + +Libgfortran is free software; you can redistribute it and/or +modify it under the terms of the GNU General Public +License as published by the Free Software Foundation; either +version 3 of the License, or (at your option) any later version. + +Libgfortran is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR a PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +Under Section 7 of GPL version 3, you are granted additional +permissions described in the GCC Runtime Library Exception, version +3.1, as published by the Free Software Foundation. + +You should have received a copy of the GNU General Public License and +a copy of the GCC Runtime Library Exception along with this program; +see the files COPYING3 and COPYING.RUNTIME respectively. If not, see +<http://www.gnu.org/licenses/>. */ + +/* This implementation of ERFC_SCALED is based on the netlib algorithm + available at http://www.netlib.org/specfun/erf */ + +#define TYPE KIND_SUFFIX(GFC_REAL_,KIND) +#define CONCAT(x,y) x ## y +#define KIND_SUFFIX(x,y) CONCAT(x,y) + +#if (KIND == 4) + +# define EXP(x) expf(x) +# define TRUNC(x) truncf(x) + +#elif (KIND == 8) + +# define EXP(x) exp(x) +# define TRUNC(x) trunc(x) + +#elif (KIND == 10) || (KIND == 16 && defined(GFC_REAL_16_IS_LONG_DOUBLE)) + +# ifdef HAVE_EXPL +# define EXP(x) expl(x) +# endif +# ifdef HAVE_TRUNCL +# define TRUNC(x) truncl(x) +# endif + +#elif (KIND == 16 && defined(GFC_REAL_16_IS_FLOAT128)) + +# define EXP(x) expq(x) +# define TRUNC(x) truncq(x) + +#else + +# error "What exactly is it that you want me to do?" + +#endif + +#if defined(EXP) && defined(TRUNC) + +extern TYPE KIND_SUFFIX(erfc_scaled_r,KIND) (TYPE); +export_proto(KIND_SUFFIX(erfc_scaled_r,KIND)); + +TYPE +KIND_SUFFIX(erfc_scaled_r,KIND) (TYPE x) +{ + /* The main computation evaluates near-minimax approximations + from "Rational Chebyshev approximations for the error function" + by W. J. Cody, Math. Comp., 1969, PP. 631-638. This + transportable program uses rational functions that theoretically + approximate erf(x) and erfc(x) to at least 18 significant + decimal digits. The accuracy achieved depends on the arithmetic + system, the compiler, the intrinsic functions, and proper + selection of the machine-dependent constants. */ + + int i; + TYPE del, res, xden, xnum, y, ysq; + +#if (KIND == 4) + static TYPE xneg = -9.382, xsmall = 5.96e-8, + xbig = 9.194, xhuge = 2.90e+3, xmax = 4.79e+37; +#else + static TYPE xneg = -26.628, xsmall = 1.11e-16, + xbig = 26.543, xhuge = 6.71e+7, xmax = 2.53e+307; +#endif + +#define SQRPI ((TYPE) 0.56418958354775628695L) +#define THRESH ((TYPE) 0.46875L) + + static TYPE a[5] = { 3.16112374387056560l, 113.864154151050156l, + 377.485237685302021l, 3209.37758913846947l, 0.185777706184603153l }; + + static TYPE b[4] = { 23.6012909523441209l, 244.024637934444173l, + 1282.61652607737228l, 2844.23683343917062l }; + + static TYPE c[9] = { 0.564188496988670089l, 8.88314979438837594l, + 66.1191906371416295l, 298.635138197400131l, 881.952221241769090l, + 1712.04761263407058l, 2051.07837782607147l, 1230.33935479799725l, + 2.15311535474403846e-8l }; + + static TYPE d[8] = { 15.7449261107098347l, 117.693950891312499l, + 537.181101862009858l, 1621.38957456669019l, 3290.79923573345963l, + 4362.61909014324716l, 3439.36767414372164l, 1230.33935480374942l }; + + static TYPE p[6] = { 0.305326634961232344l, 0.360344899949804439l, + 0.125781726111229246l, 0.0160837851487422766l, + 0.000658749161529837803l, 0.0163153871373020978l }; + + static TYPE q[5] = { 2.56852019228982242l, 1.87295284992346047l, + 0.527905102951428412l, 0.0605183413124413191l, + 0.00233520497626869185l }; + + y = (x > 0 ? x : -x); + if (y <= THRESH) + { + ysq = 0; + if (y > xsmall) + ysq = y * y; + xnum = a[4]*ysq; + xden = ysq; + for (i = 0; i <= 2; i++) + { + xnum = (xnum + a[i]) * ysq; + xden = (xden + b[i]) * ysq; + } + res = x * (xnum + a[3]) / (xden + b[3]); + res = 1 - res; + res = EXP(ysq) * res; + return res; + } + else if (y <= 4) + { + xnum = c[8]*y; + xden = y; + for (i = 0; i <= 6; i++) + { + xnum = (xnum + c[i]) * y; + xden = (xden + d[i]) * y; + } + res = (xnum + c[7]) / (xden + d[7]); + } + else + { + res = 0; + if (y >= xbig) + { + if (y >= xmax) + goto finish; + if (y >= xhuge) + { + res = SQRPI / y; + goto finish; + } + } + ysq = ((TYPE) 1) / (y * y); + xnum = p[5]*ysq; + xden = ysq; + for (i = 0; i <= 3; i++) + { + xnum = (xnum + p[i]) * ysq; + xden = (xden + q[i]) * ysq; + } + res = ysq *(xnum + p[4]) / (xden + q[4]); + res = (SQRPI - res) / y; + } + +finish: + if (x < 0) + { + if (x < xneg) + res = __builtin_inf (); + else + { + ysq = TRUNC (x*((TYPE) 16))/((TYPE) 16); + del = (x-ysq)*(x+ysq); + y = EXP(ysq*ysq) * EXP(del); + res = (y+y) - res; + } + } + return res; +} + +#endif + +#undef EXP +#undef TRUNC + +#undef CONCAT +#undef TYPE +#undef KIND_SUFFIX |