[Bug fortran/65757] New: gfortran gives incorrect result for anint with real*16 argument

[Bug fortran/65757] New: gfortran gives incorrect result for anint with real*16 argument

[bugs at dhbailey dot com]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757

            Bug ID: 65757
           Summary: gfortran gives incorrect result for anint with real*16
                    argument
           Product: gcc
           Version: 4.9.0
            Status: UNCONFIRMED
          Severity: major
          Priority: P3
         Component: fortran
          Assignee: unassigned at gcc dot gnu.org
          Reporter: bugs at dhbailey dot com

gfortran gives an incorrect result for the anint (nearest whole number)
function for certain real*16 arguments.  Here is a simple example:

program anintbug
implicit none
real (kind (0.q0)) q1, q2, q3

q1 = 233181505644407.99996948242187500000q0
q2 = anint (q1)
q3 = 233181505644408.00000000000000000000q0

write (6, '(f50.25)') q1, q2, q3
write (6, '(z35)') q1, q2, q3
stop
end

When compiled with "gfortran anintbug.f90", this program produces the following
output.  The second result should equal the third result, but does not.

         233181505644407.9999694824218750000000000
         233181505644407.9999694824218750000000000
         233181505644408.0000000000000000000000000
   402EA827999FCEEFFFFC000000000000
   402EA827999FCEEFFFFC000000000000
   402EA827999FCEF00000000000000000

[Bug libquadmath/65757] gfortran gives incorrect result for anint with real*16 argument

[kargl at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757

kargl at gcc dot gnu.org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |kargl at gcc dot gnu.org
          Component|fortran                     |libquadmath
      Known to fail|                            |4.6.4, 4.7.4, 4.8.5, 4.9.3,
                   |                            |5.0, 6.0
           Severity|major                       |normal

--- Comment #1 from kargl at gcc dot gnu.org ---
Using -fdump-tree-original (and removing the IO code), 
one gets

anintbug ()
{
  real(kind=16) q1;
  real(kind=16) q2;
  real(kind=16) q3;

  q1 = 2.33181505644407999969482421875e+14;
  q2 = roundq (q1);
  q3 = 2.33181505644408e+14;

}

So, gfortran is calling the libquadmath routine roundq.
I've changed the component from fortran to libquadmath.
This also appears to fail with all versions of gfortran.

[Bug libquadmath/65757] gfortran gives incorrect result for anint with real*16 argument

[sgk at troutmask dot apl.washington.edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757

--- Comment #2 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Tue, Apr 14, 2015 at 01:38:00AM +0000, kargl at gcc dot gnu.org wrote:
> Using -fdump-tree-original (and removing the IO code), 
> one gets
> 
> anintbug ()
> {
>   real(kind=16) q1;
>   real(kind=16) q2;
>   real(kind=16) q3;
> 
>   q1 = 2.33181505644407999969482421875e+14;
>   q2 = roundq (q1);
>   q3 = 2.33181505644408e+14;
> 
> }
> 
> So, gfortran is calling the libquadmath routine roundq.
> I've changed the component from fortran to libquadmath.
> This also appears to fail with all versions of gfortran.
> 

The following patch appears to fix this issue, but I
do not know the internals of __float128_t so the bit
twiddling may not be right.

Index: libquadmath/math/roundq.c
===================================================================
--- libquadmath/math/roundq.c    (revision 222054)
+++ libquadmath/math/roundq.c    (working copy)
@@ -69,7 +69,7 @@ roundq (__float128 x)
     }
   else
     {
-      uint64_t i = -1ULL >> (j0 - 48);
+      uint64_t i = -1ULL >> (j0 - 47);
       if ((i1 & i) == 0)
     /* X is integral.  */
     return x;
@@ -77,7 +77,7 @@ roundq (__float128 x)
       if (huge + x > 0.0)
     {
       /* Raise inexact if x != 0.  */
-      uint64_t j = i1 + (1LL << (111 - j0));
+      uint64_t j = i1 + (1LL << (110 - j0));
       if (j < i1)
         i0 += 1;
       i1 = j;

[Bug libquadmath/65757] gfortran gives incorrect result for anint with real*16 argument

[bugs at dhbailey dot com]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757

--- Comment #3 from dhbbugs <bugs at dhbailey dot com> ---
Has anyone else been able to confirm that Steve Kargl's fix will work?  What
has to happen to move this fix into the production code?

[Bug libquadmath/65757] gfortran gives incorrect result for anint with real*16 argument

[jvdelisle at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757

Jerry DeLisle <jvdelisle at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |jvdelisle at gcc dot gnu.org

--- Comment #4 from Jerry DeLisle <jvdelisle at gcc dot gnu.org> ---
Typically you need a libquadmath maintainer to review the patch and approve it.

I did try the patch with the small test case here, but it probably needs to get
tested for other cases.

Since it is a fairly serious wrong code bug, it would first need to go into GCC
trunk (6.0) and then after some proving time, it could be back ported. The
patch seems simple enough that one could go back a few versions I suppose.

[Bug libquadmath/65757] gfortran gives incorrect result for anint with real*16 argument

[sgk at troutmask dot apl.washington.edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757

--- Comment #5 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Wed, Apr 22, 2015 at 11:19:13PM +0000, bugs at dhbailey dot com wrote:
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757
> 
> --- Comment #3 from dhbbugs <bugs at dhbailey dot com> ---
> Has anyone else been able to confirm that Steve Kargl's fix will work?  What
> has to happen to move this fix into the production code?
> 

I haven't submitted the patch to gcc-patches@, so I suspect
no one has looked at it (especially with gcc 5.1 being 
released).  I simply attached the patch to the audit trail
of the bug report because I haven't been able to convince myself
that it is correct.

[Bug libquadmath/65757] gfortran gives incorrect result for anint with real*16 argument

[jvdelisle at gcc dot gnu.org]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757

--- Comment #6 from Jerry DeLisle <jvdelisle at gcc dot gnu.org> ---
I am thinking its not working, but maybe I do not understand the function.  I
created a loop like this (with the patch):

q1 = 233181505644407.99999999999999999999q0
q2 = anint (q1)
q3 = 233181505644408.00000000000000000000q0
do i=1,1000
  q1=q1+.02400039q0
  q2=anint(q1)
  write (6, '(f50.25)') q1, q2
  write (6, '(z35)') q1, q2
end do

Scrolling through the results, I spot this:

         233181505644431.7363857099999999995104932
         233181505644431.5000000000000000000000000
   402EA827999FCF1F79078C3B622BC417
   402EA827999FCF1F0000000000000000

[Bug libquadmath/65757] gfortran gives incorrect result for anint with real*16 argument

[dominiq at lps dot ens.fr]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757

Dominique d'Humieres <dominiq at lps dot ens.fr> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |burnus at gcc dot gnu.org,
                   |                            |jakub at gcc dot gnu.org

--- Comment #7 from Dominique d'Humieres <dominiq at lps dot ens.fr> ---
CCing libquadmath maintainers.

[Bug libquadmath/65757] gfortran gives incorrect result for anint with real*16 argument

[bugs at dhbailey dot com]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757

--- Comment #8 from dhbbugs <bugs at dhbailey dot com> ---
Jerry DeLisl'e output is certainly not correct -- anint should invariably
return the nearest whole number.  It should be the equivalent of this code:

if (x >= 0.0) then
  anint = aint (x + 0.5)
else
  anint = aint (x - 0.5)
endif

Here aint is the truncated whole number function (a Fortran intrinsic).

[Bug libquadmath/65757] gfortran gives incorrect result for anint with real*16 argument

[joseph at codesourcery dot com]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757

--- Comment #9 from joseph at codesourcery dot com <joseph at codesourcery dot com> ---
Fixed in glibc (commit 7d0b2575416aec2717e8665287d0ab77826a0ade).  I'd 
advise merging to trunk GCC libquadmath all relevant glibc changes since 
the last merges in 2012, rather than cherry-picking just this fix 
(although if you wish to fix things on release branches, cherry-picking 
just the critical fixes would be appropriate there).

[Bug libquadmath/65757] gfortran gives incorrect result for anint with real*16 argument

[sgk at troutmask dot apl.washington.edu]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757

--- Comment #10 from Steve Kargl <sgk at troutmask dot apl.washington.edu> ---
On Tue, Apr 28, 2015 at 05:32:11PM +0000, joseph at codesourcery dot com wrote:
> https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757
> 
> --- Comment #9 from joseph at codesourcery dot com <joseph at codesourcery dot com> ---
> Fixed in glibc (commit 7d0b2575416aec2717e8665287d0ab77826a0ade).  I'd 
> advise merging to trunk GCC libquadmath all relevant glibc changes since 
> the last merges in 2012, rather than cherry-picking just this fix 
> (although if you wish to fix things on release branches, cherry-picking 
> just the critical fixes would be appropriate there).
> 

Short of grabbing the entire glibc repository, how does
one grab the current libquadmath code?  In particular, 
gitweb suggests that there isn't a component of glibc
with the name libquadmath.

[Bug libquadmath/65757] gfortran gives incorrect result for anint with real*16 argument

[joseph at codesourcery dot com]

https://gcc.gnu.org/bugzilla/show_bug.cgi?id=65757

--- Comment #11 from joseph at codesourcery dot com <joseph at codesourcery dot com> ---
I don't know what process Jakub and Tobias used (a) to identify relevant 
files / changes in glibc and (b) to make all the changes to operate on 
__float128 rather than long double.  Roughly, it's 
sysdeps/ieee754/ldbl-128, *plus* strtod-related files from stdlib/, *plus* 
printf-related files from stdio-common, *plus* some functions (especially 
<complex.h> functions) from math/, that are used in libquadmath and so may 
need updating for changes made in glibc.