PR fortran/33197

gcc/fortran/
	* intrinsic.c (add_functions): Modify intrinsics ACOSH, ASINH,
	ATANH, ERF, ERFC and GAMMA. Add intrinsics BESSEL_{J,Y}{0,1,N},
	ERFC_SCALED, LOG_GAMMA and HYPOT.
	* intrinsic.h (gfc_check_hypot, gfc_simplify_hypot,
	gfc_resolve_hypot): New prototypes.
	* mathbuiltins.def: Add HYPOT builtin. Make complex versions of
	ACOSH, ASINH and ATANH available.
	* gfortran.h (GFC_ISYM_ERFC_SCALED, GFC_ISYM_HYPOT): New values.
	* lang.opt: Add -std=f2008 option.
	* libgfortran.h: Define GFC_STD_F2008.
	* lang-specs.h: Add .f08 and .F08 file suffixes.
	* iresolve.c (gfc_resolve_hypot): New function.
	* parse.c (parse_contained): Allow empty CONTAINS for Fortran 2008.
	* check.c (gfc_check_hypot): New function.
	* trans-intrinsic.c (gfc_intrinsic_map): Define ERFC_SCALE builtin.
	* options.c (set_default_std_flags): Allow Fortran 2008 by default.
	(form_from_filename): Add .f08 suffix.
	(gfc_handle_option): Handle -std=f2008 option.
	* simplify.c (gfc_simplify_hypot): New function.
	* gfortran.texi: Document Fortran 2008 status and file extensions.
	* intrinsic.texi: Document new BESSEL_{J,Y}{0,1,N} intrinsics,
	as well as HYPOT and ERFC_SCALED. Update documentation of ERF,
	ERFC, GAMMA, LGAMMA, ASINH, ACOSH and ATANH.
	* invoke.texi: Document the new -std=f2008 option.

libgomp/
	* testsuite/libgomp.fortran/fortran.exp: Add .f08 and
	.F08 file suffixes.

gcc/testsuite/
	* gfortran.dg/gomp/gomp.exp: Add .f08 and .F08 file suffixes.
	* gfortran.dg/dg.exp: Likewise.
	* gfortran.dg/vect/vect.exp: Likewise.
	* gfortran.fortran-torture/execute/execute.exp: Likewise.
	* gfortran.fortran-torture/compile/compile.exp: Likewise.
	* gfortran.dg/gamma_1.f90: Also check log_gamma.
	* gfortran.dg/invalid_contains_1.f90: Remove warning about
	empty CONTAINS.
	* gfortran.dg/gamma_2.f90: Add a few error messages.
	* gfortran.dg/invalid_contains_2.f90: Remove warning about
	empty CONTAINS.
	* gfortran.dg/gamma_3.f90: Adjust error message.
	* gfortran.dg/gamma_4.f90: Test for log_gamma instead of lgamma.
	* gfortran.dg/bind_c_usage_9.f03: Adjust error messages.
	* gfortran.dg/bessel_1.f90: New test.
	* gfortran.dg/recursive_check_3.f90: Remove warnings.
	* gfortran.dg/besxy.f90: Also check for new F2008 intrinsics.
	* gfortran.dg/derived_function_interface_1.f90: Remove warning.
	* gfortran.dg/contains_empty_1.f03: New test.
	* gfortran.dg/erfc_scaled_1.f90: New test.
	* gfortran.dg/hypot_1.f90: New test.
	* gfortran.dg/contains_empty_2.f03: New test.

libgfortran/
	* intrinsics/erfc_scaled_inc.c: New file.
	* intrinsics/erfc_scaled.c: New file.
	* gfortran.map (GFORTRAN_1.0): Add _gfortran_erfc_scaled_r*.
	* Makefile.am: Add intrinsics/erfc_scaled.c.
	* config.h.in: Regenerate.
	* configure: Regenerate.
	* Makefile.in: Regenerate.


git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@132846 138bc75d-0d04-0410-961f-82ee72b054a4

1 files changed, 175 insertions, 0 deletions

diff --git a/libgfortran/intrinsics/erfc_scaled_inc.c b/libgfortran/intrinsics/erfc_scaled_inc.c
new file mode 100644
index 00000000000..fab815584ca
--- /dev/null
+++ b/libgfortran/intrinsics/erfc_scaled_inc.c
@@ -0,0 +1,175 @@
+/* Implementation of the ERFC_SCALED intrinsic, to be included by erfc_scaled.c
+   Copyright (c) 2008 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 2 of the License, or (at your option) any later version.
+
+In addition to the permissions in the GNU General Public License, the
+Free Software Foundation gives you unlimited permission to link the
+compiled version of this file into combinations with other programs,
+and to distribute those combinations without any restriction coming
+from the use of this file.  (The General Public License restrictions
+do apply in other respects; for example, they cover modification of
+the file, and distribution when not linked into a combine
+executable.)
+
+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.
+
+You should have received a copy of the GNU General Public
+License along with libgfortran; see the file COPYING.  If not,
+write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+Boston, MA 02110-1301, USA.  */
+
+/* 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)
+#else
+# define EXP(x) expl(x)
+# define TRUNC(x) truncl(x)
+#endif
+
+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;
+}
+
+#undef EXP
+#undef TRUNC
+
+#undef CONCAT
+#undef TYPE
+#undef KIND_SUFFIX