re PR libfortran/49024 (REAL*16 ERFC_SCALED inaccuracy)

	PR libfortran/49024

	* intrinsics/erfc_scaled.c (erfc_scaled_r16): New function.
	* intrinsics/erfc_scaled_inc.c: Do not provide quadruple
	precision variant.

	* gfortran.dg/erf_3.F90: New file.

From-SVN: r205151

diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog
index d47ec92a4be2f2c1f966cd3ac25640677ad4fc47..694fb831fdaa8a1f98f60121c27fc5f8e4ff89d8 100644 (file)
--- a/gcc/testsuite/ChangeLog
+++ b/gcc/testsuite/ChangeLog
@@ -1,3 +1,8 @@
+2013-11-20  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
+
+       PR libfortran/49024
+       * gfortran.dg/erf_3.F90: New file.
+
 2013-11-20  Bill Schmidt  <wschmidt@linux.vnet.ibm.com>
 
        * gcc.target/powerpc/pr48258-1.c: Skip for little endian.

diff --git a/gcc/testsuite/gfortran.dg/erf_3.F90 b/gcc/testsuite/gfortran.dg/erf_3.F90
new file mode 100644 (file)
index 0000000..fcff41f
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/erf_3.F90
@@ -0,0 +1,44 @@
+! { dg-do run }
+! { dg-options "-fno-range-check -ffree-line-length-none -O0" }
+! { dg-add-options ieee }
+!
+! Check that simplification functions and runtime library agree on ERF,
+! ERFC and ERFC_SCALED, for quadruple-precision.
+
+program test
+  implicit none
+
+  real(kind=16) :: x16
+
+#define CHECK(a) \
+  x16 = a ; \
+  call check(erf(real(a,kind=16)), erf(x16)) ; \
+  call check(erfc(real(a,kind=16)), erfc(x16)) ; \
+  call check(erfc_scaled(real(a,kind=16)), erfc_scaled(x16))
+
+  CHECK(0.0)
+  CHECK(0.9)
+  CHECK(1.9)
+  CHECK(10.)
+  CHECK(11.)
+  CHECK(12.)
+  CHECK(13.)
+  CHECK(14.)
+  CHECK(49.)
+  CHECK(190.)
+
+  CHECK(-0.0)
+  CHECK(-0.9)
+  CHECK(-1.9)
+  CHECK(-19.)
+  CHECK(-190.)
+
+contains
+
+  subroutine check (a, b)
+    real(kind=16), intent(in) :: a, b
+    print *, abs(a-b) / spacing(a)
+    if (abs(a - b) > 10 * spacing(a)) call abort
+  end subroutine
+
+end program test

diff --git a/libgfortran/ChangeLog b/libgfortran/ChangeLog
index fcbc54834722911e29ba571307f1ef37e55d5f0d..eeeaa04b3c3ca389e8ce39c7b5476f76ba444a6d 100644 (file)
--- a/libgfortran/ChangeLog
+++ b/libgfortran/ChangeLog
@@ -1,3 +1,10 @@
+2013-11-20  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
+
+       PR libfortran/49024
+       * intrinsics/erfc_scaled.c (erfc_scaled_r16): New function.
+       * intrinsics/erfc_scaled_inc.c: Do not provide quadruple
+       precision variant.
+
 2013-11-18  Francois-Xavier Coudert  <fxcoudert@gcc.gnu.org>
 
        PR libfortran/51828

diff --git a/libgfortran/intrinsics/erfc_scaled.c b/libgfortran/intrinsics/erfc_scaled.c
index 1c58a08938e0d95108d32a3bc7d8c5b7148c9469..1f8c778eb651b638f7b1a3b5614966537baca19c 100644 (file)
--- a/libgfortran/intrinsics/erfc_scaled.c
+++ b/libgfortran/intrinsics/erfc_scaled.c
@@ -45,8 +45,59 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
 #include "erfc_scaled_inc.c"
 #endif
 
-#ifdef HAVE_GFC_REAL_16
+#if defined(HAVE_GFC_REAL_16) && defined(GFC_REAL_16_IS_LONG_DOUBLE)
 #undef KIND
 #define KIND 16
 #include "erfc_scaled_inc.c"
 #endif
+
+
+/* For quadruple-precision (__float128), netlib's implementation is
+   not accurate enough.  We provide another one.  */
+
+
+extern GFC_REAL_16 erfc_scaled_r16 (GFC_REAL_16);
+export_proto(erfc_scaled_r16);
+
+
+GFC_REAL_16
+erfc_scaled_r16 (GFC_REAL_16 x)
+{
+  if (x < -106.566990228185312813205074546585730Q)
+    {
+      return __builtin_infq();
+    }
+  if (x < 12)
+    {
+      /* Compute directly as ERFC_SCALED(x) = ERFC(x) * EXP(X**2).
+        This is not perfect, but much better than netlib.  */
+      return erfcq(x) * expq(x * x);
+    }
+  else
+    {
+      /* Calculate ERFC_SCALED(x) using a power series in 1/x:
+        ERFC_SCALED(x) = 1 / (x * sqrt(pi))
+                        * (1 + Sum_n (-1)**n * (1 * 3 * 5 * ... * (2n-1))
+                                             / (2 * x**2)**n)
+       */
+      GFC_REAL_16 sum = 0, oldsum;
+      GFC_REAL_16 inv2x2 = 1 / (2 * x * x);
+      GFC_REAL_16 fac = 1;
+      int n = 1;
+
+      while (n < 200)
+       {
+         fac *= - (2*n - 1) * inv2x2;
+         oldsum = sum;
+         sum += fac;
+
+         if (sum == oldsum)
+           break;
+
+         n++;
+       }
+
+      return (1 + sum) / x * (M_2_SQRTPIq / 2);
+    }
+}
+

diff --git a/libgfortran/intrinsics/erfc_scaled_inc.c b/libgfortran/intrinsics/erfc_scaled_inc.c
index 57a6b71f995dd85bb3bf29e14c34cfe199fafd52..107d91a6c9db2c2c03cfcab637901ba349bcf3ae 100644 (file)
--- a/libgfortran/intrinsics/erfc_scaled_inc.c
+++ b/libgfortran/intrinsics/erfc_scaled_inc.c
@@ -48,11 +48,6 @@ see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
 #  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?"