Author: fxcoudert
Date: Thu Nov 21 08:45:00 2013
New Revision: 205193

URL: http://gcc.gnu.org/viewcvs?rev=205193&root=gcc&view=rev
Log:
	PR libfortran/59227
	* intrinsics/erfc_scaled.c (erfc_scaled_r16): Don't define if
	__float128 is not available.

Modified:
    trunk/libgfortran/ChangeLog
    trunk/libgfortran/intrinsics/erfc_scaled.c

Fixed on trunk. Sorry, and thanks Andreas for the report!

That's a different bug.

What does this output, when compiled with "-fno-range-check -ffree-line-length-none -O0"?

$ cat test.f90
program test
  real(kind=16) :: x
  x = 12
  print *, erfc_scaled(real(12,kind=16))
  print *, erfc_scaled(x)
end program test

   4.68542210148937626195884133993966578E-0002
  -1.87533922948603221391158058278771595E+0920

(In reply to Andreas Schwab from comment #5)
>    4.68542210148937626195884133993966578E-0002
>   -1.87533922948603221391158058278771595E+0920

And this (at -O2)?

__float128 foo (__float128 x)
{
  __float128 sum = 0, oldsum;
  __float128 inv2x2 = 1 / (2 * x * x);
  __float128 fac = 1;
  int n = 1;

  while (n < 200)
  {
    fac *= - (2*n - 1) * inv2x2;
    oldsum = sum;
    sum += fac;

    __builtin_printf ("%lg\n", (double) fac);
    if (sum == oldsum)
      break;

    n++;
  }

  return (1 + sum) / x * (1.1283791670955125738961589031215452Q / 2);
}

int main (void)
{
  __float128 x = 12;
  x = foo(x);
  __builtin_printf ("%lg\n", (double) x);
}

-0.00347222
3.6169e-05
-6.27934e-07
1.52623e-08
-4.76946e-10
1.82167e-11
-8.22281e-13
4.28272e-14
-2.52799e-15
1.66777e-16
-1.21608e-17
9.71178e-19
-8.43037e-20
7.90347e-21
-7.95835e-22
8.56628e-23
-9.81553e-24
1.19286e-24
-1.53249e-25
2.07525e-26
-2.95435e-27
4.41101e-28
-6.8922e-29
1.12477e-29
-1.91367e-30
3.38879e-31
-6.23632e-32
1.19096e-32
-2.35711e-33
4.82881e-34
-1.02277e-34
2.23731e-35
-5.04948e-36
1.17471e-36
-2.8144e-37
6.93827e-38
0.0468542

So the calculation inside __gfortran_erfc_scaled_r16 is correct, but the result is printed incorrectly. Can you try this:

program test
  real(kind=16) :: x
  x = 10.9_16 ; print *, erfc_scaled(x)
  x = 11.9_16 ; print *, erfc_scaled(x)
  x = 12.0_16 ; print *, erfc_scaled(x)
  x = 12.1_16 ; print *, erfc_scaled(x)
  x = 13.1_16 ; print *, erfc_scaled(x)
  x = 14.1_16 ; print *, erfc_scaled(x)
end program test

Also, can you confirm that binary128 is not the same as long double on ia64?

(I've sent a request to reopen my old account on the compile farm… If this drags out, I'll investigate on my own there. But for now, I don't have access to an ia64.)

Looks like inv2x2 is wrong in _gfortran_erfc_scaled_r16:

(gdb) i reg r42 r43
r42            0x0      0
r43            0x4007200000000000       4613691527636451328

   5.15453772226367323012693682499968245E-0002
   4.72452324840876699523132987063922456E-0002
  -1.87533922948603221391158058278771595E+0920
  -5.09942677661379357690770320803840300E+0921
  -2.70969365071138516780674252318715535E+0935
  -1.40653905502094549174911373980686654E+0948

So, it looks like the same code that works fine in an external C program, is miscompiled in libgfortran's _gfortran_erfc_scaled_r16. Do you agree?

Can you remove the __builtin_printf call inside foo() in comment #6, and compile with the same flags as used to compile libgfortran? For me, that would be "-std=gnu99 -fcx-fortran-rules -ffunction-sections -fdata-sections -g -O2" Just to be sure…

1 / (2 * x * x) is calling __divtf3 from glibc, which operates on long double, not quad float.  The function from libgcc.a should have been called instead.

(In reply to Andreas Schwab from comment #12)
> 1 / (2 * x * x) is calling __divtf3 from glibc, which operates on long
> double, not quad float.  The function from libgcc.a should have been called
> instead.

Because ia64 uses TFmode for both long double and __float128, while i386 uses XFmode for long double, is that it?

I'm not sure what I can do, then… and don't understand why this does not happen in the standalone test case in comment #6.

I think this can be considered a glibc bug, which should not have defined __divtf3 in the first place (the function should have been called __divxf3), and  the test should be XFAILd on ia64.

(In reply to Andreas Schwab from comment #14)
> I think this can be considered a glibc bug, which should not have defined
> __divtf3 in the first place (the function should have been called __divxf3),
> and  the test should be XFAILd on ia64.

Is "ia64-*-linux" the right triplet? If you can confirm that the following patch works as expect on your system, I'll commit it:

Index: erf_3.F90
===================================================================
--- erf_3.F90	(revision 205151)
+++ erf_3.F90	(working copy)
@@ -1,20 +1,29 @@
-! { dg-do run }
+! { dg-do run { xfail spu-*-* ia64-*-linux } }
 ! { 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.
+!
+! XFAILed for SPU targets because our library implementation of
+! the double-precision erf/erfc functions is not accurate enough.
+!
+! XFAILed for IA64 Linux because of a glibc bug:
+! http://gcc.gnu.org/bugzilla/show_bug.cgi?id=59227
 
 program test
+  use, intrinsic :: iso_fortran_env
   implicit none
 
-  real(kind=16) :: x16
+  ! QP will be the largest supported real kind, possibly real(kind=16)
+  integer, parameter :: qp = real_kinds(ubound(real_kinds,dim=1))
+  real(kind=qp) :: x
 
 #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))
+  x = a ; \
+  call check(erf(real(a,kind=qp)), erf(x)) ; \
+  call check(erfc(real(a,kind=qp)), erfc(x)) ; \
+  call check(erfc_scaled(real(a,kind=qp)), erfc_scaled(x))
 
   CHECK(0.0)
   CHECK(0.9)
@@ -36,7 +45,7 @@ program test
 contains
 
   subroutine check (a, b)
-    real(kind=16), intent(in) :: a, b
+    real(kind=qp), intent(in) :: a, b
     print *, abs(a-b) / spacing(a)
     if (abs(a - b) > 10 * spacing(a)) call abort
   end subroutine

Perhaps a suitable workaround would be to force linking the right __divtf3 into libgfortran.

ia64 was using TFmode for 80 bit long double before r72916, so the __divtf3 name sticked for compatibility.

Author: fxcoudert
Date: Thu Nov 21 11:37:07 2013
New Revision: 205210

URL: http://gcc.gnu.org/viewcvs?rev=205210&root=gcc&view=rev
Log:
	PR libfortran/59227
	* gfortran.dg/erf_3.F90: XFAIL on spu-* and ia64-*-linux*.
	Make more generic for other platforms.

Modified:
    trunk/gcc/testsuite/ChangeLog
    trunk/gcc/testsuite/gfortran.dg/erf_3.F90

(In reply to Andreas Schwab from comment #16)
> Perhaps a suitable workaround would be to force linking the right __divtf3
> into libgfortran.

I'm not sure how to force linking in libgfortran, but if it helps reliably fix quad-precision math on ia64, that sounds like a good option.

Meanwhile, I've xfail'ed the test (also on spu-* due to a library issue there, similar to erf_2.F90).

It's not really a glibc bug, but a compatilibity wart.

> I'm not sure how to force linking in libgfortran

Add libgcc/divtf3_s.o to the link line.

Actually it appears to be a libgcc bug.

I saw in your daily tests that the XFAIL indeed silences the failure. Closing as fixed, will un-XFAIL when PR59230 is fixed.

The master branch has been updated by Andreas Schwab <schwab@gcc.gnu.org>:

https://gcc.gnu.org/g:0ecf5229723ec99e6b5099dd68d48bd925da6b0d

commit r11-895-g0ecf5229723ec99e6b5099dd68d48bd925da6b0d
Author: Andreas Schwab <schwab@suse.de>
Date:   Thu Nov 21 15:15:40 2013 +0100

    Missing __divtf3@@GCC_4.4.0 on ia64
    
    gcc/testsuite/
            PR libfortran/59227
            * gfortran.dg/erf_3.F90: Remove XFAIL on ia64-*-linux*.
    
    libgcc/
            PR target/59230
            PR libfortran/59227
            * config/ia64/t-softfp-compat (softfp_file_list): Filter out
            soft-fp/divtf3.c.
            (LIB2ADD): Add config/ia64/divtf3.c.
            * config/ia64/divtf3.c: New file.