[Bug math/4997] New: lround(nexafter(0.5,-1)) has incorrect value

[Bug math/4997] New: lround(nexafter(0.5,-1)) has incorrect value

[fxcoudert at gcc dot gnu dot org]

On powerpc-unknown-linux-gnu, with GNU C Library development release version 2.4
(20060617), lround(nextafter(0.5,-1)) returns 1 instead of the expected 0.


PS: it appears that powerpc glibc-2.6 is calculating lround(x) by converting
(x+0.5) into an integer (thanks to David Edelsohn for pointing that to me):

    fcmpu   cr6,fp1,fp12    /* if (x > 0.0)  */
    ble-    cr6,.L4
    fadd    fp1,fp1,fp10    /* x+= 0.5;  */
    fctiwz  fp2,fp1         /* Convert To Integer DW lround toward 0.  */
.L4:
    fsub    fp1,fp1,fp10    /* x-= 0.5;  */

I think 0.5 needs to be changed into nextafter(0.5,-1). Details of why can be
found here: http://gcc.gnu.org/ml/fortran/2005-04/msg00139.html

-- 
           Summary: lround(nexafter(0.5,-1)) has incorrect value
           Product: glibc
           Version: 2.4
            Status: NEW
          Severity: normal
          Priority: P2
         Component: math
        AssignedTo: aj at suse dot de
        ReportedBy: fxcoudert at gcc dot gnu dot org
                CC: glibc-bugs at sources dot redhat dot com
 GCC build triplet: powerpc-unknown-linux-gnu
  GCC host triplet: powerpc-unknown-linux-gnu
GCC target triplet: powerpc-unknown-linux-gnu


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[fxcoudert at gcc dot gnu dot org]

------- Additional Comments From fxcoudert at gcc dot gnu dot org  2007-09-14 12:36 -------
(In reply to comment #0)
> On powerpc-unknown-linux-gnu

Some more testing indicates that powerpc64-unknown-linux-gnu is not affected by
the problem.

-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[fxcoudert at gcc dot gnu dot org]

------- Additional Comments From fxcoudert at gcc dot gnu dot org  2007-09-14 13:54 -------
Yet some more information:

  -- ppc64-linux / glibc-2.4-20070330 works OK, with both -m32 and -m64
  -- ppc-linux / glibc-2.6.1 fails, with both -m32 and -m64
  -- ppc-linux / glibc-2.4-20060617 fails with -m32 (don't know about -m64)

-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[rsa at us dot ibm dot com]

------- Additional Comments From rsa at us dot ibm dot com  2007-10-22 22:52 -------
We've identified a bug in the powerpc version of lround whereby if x == 0.5 or
-0.5 exactly the fadd in lround will round improperly due to the guard bits
forcing a shift and round.

Likewise there is a similar case in llround whereby if x is sufficiently large
(max double or min double) and an fadd of 0.5 or -0.5 is performed in llround
the guard bits force a shift and the resultant value will be incorrect.

I'll work on a fix for both cases.

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
         AssignedTo|aj at suse dot de           |rsa at us dot ibm dot com
             Status|NEW                         |ASSIGNED


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[rsa at us dot ibm dot com]

------- Additional Comments From rsa at us dot ibm dot com  2007-11-01 22:13 -------
Created an attachment (id=2072)
 --> (http://sourceware.org/bugzilla/attachment.cgi?id=2072&action=view)
libm-test.inc patch to identify lround and llround erroneous round.

This is a patch to libm-test.inc which will identify the erroneous lround and
llround cases mentioned in this bug description.  We'll identify a fix and
include these test-cases patch in the eventually fix.

-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[rsa at us dot ibm dot com]

------- Additional Comments From rsa at us dot ibm dot com  2007-11-08 18:53 -------
The following asm stub resolves the issue by returning 'zero' if 0.5 > x > -0.5.  

        fsub    fp10,fp0,fp0    /* Compute 0.0.  */
        fabs    fp2, fp1        /* Get the absolute value of x.  */
        fcmpu   cr3, fp1, fp10  /* orig x is negative? x < 0.0  */
        fcmpu   cr6, fp2, fp0   /* if |x| < 0.5  */
        ble-    cr6,.Lretzero
        fadd    fp3,fp2,fp0     /* |x|+=0.5 bias to prepare to round.  */
        bge     cr3,.Lconvert   /* orig x was positive so don't negate ans.  */
        fnabs   fp2,fp3         /* -(|x|+=0.5)  */
.Lconvert:
        fctiwz  fp4,fp2         /* Convert to Integer DW lround toward 0.  */
        stfd    fp4,8(r1)
        nop
        nop
        nop
        lwz     r3,12(r1)       /* Load return as integer.  */
.Lout:
        addi    r1,r1,16
        blr
.Lretzero:              /* 0.5 > x > -0.5  */
        li      r3,0            /* return 0.  */
        b       .Lout

I need to identify all the places where the erroneous rounding is happening,
apply this fix and also apply a check where x > 2^52 for llround.  This may also
be an issue for round().  Patch forthcoming.


-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[fxcoudert at gcc dot gnu dot org]

------- Additional Comments From fxcoudert at gcc dot gnu dot org  2007-11-08 18:55 -------
Doesn't the same problem happen for 1.5 and 2.5?

-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[rsa at us dot ibm dot com]

------- Additional Comments From rsa at us dot ibm dot com  2007-11-09 13:45 -------
Francois-Xavier, The bug doesn't manifest for values nextafter(1.5,-1),
nextafter(2.5,-1), etc because nextafter(0.5,-1) just happens to be encoded with
all 1's in the mantissa, e.g. 0x1.fffffffffffffp-2 and the shift and carry
results in an erroneous round-up due to the fadd instruction.  All 1's in the
mantissa is not the resultant encoding for nextafter(1.5,-1)
[0x1.7ffffffffffffp+0] or others.

The other cases where there are all 1's in the mantissa are:
0x1.fffffffffffffp-1 == .999999 which should rightly round up.
0x1.fffffffffffffp+0 == 1.99999 which should rightly round up, etc.

So nextafter(+0.5,-1) is the only case for lround where all 1's in the mantissa
should round down instead of up.

We fixed a couple bugs where a stale register was being passed to fctiwz.  Also
the .Lretzero branch should be 'less than' per the comments and not 'less than
or equal' as in the previous asm code.  Finally we reordered the initial
instructions for data dependence streamlining.

        lis 9,.LC1@ha
        lfs 10,.LC1@l(9)        /* Load constant 0.5 into fpr10.  */
        fabs    2, 1            /* Get the absolute value of x.  */
        fsub    12,10,10        /* Compute 0.0.  */
        fcmpu   6, 2, 10        /* if |x| < 0.5  */
        fcmpu   3, 1, 12        /* x is negative? x < 0.0  */
        blt-    6,.Lretzero
        fadd    3,2,10          /* |x|+=0.5 bias to prepare to round.  */
        bge     3,.Lconvert     /* x is positive so don't negate x.  */
        fnabs   3,3             /* -(|x|+=0.5)  */
.Lconvert:
        fctiwz  4,3             /* Convert to Integer word lround toward 0.  */
        stfd    4,8(1)
        nop
        nop
        nop
        lwz     3,12(1)         /* Load return as integer.  */
.Lout:
        addi    1,1,16
        blr
.Lretzero:                      /* 0.5 > x > -0.5  */
        li      3,0             /* return 0.  */

I still have to add the case where x is between 2^52 and 2^53-1 to cover the
llround case.

Expect  0: Got  0 = lround(nextafter(0.5,-1) [0x1.fffffffffffffp-2]
{0.4999999999999999})
Expect  0: Got  0 = lround(nextafter(-0.5,1) [-0x1.fffffffffffffp-2]
{-0.4999999999999999})
Expect  1: Got  1 = lround(nextafter(1.5,-1) [0x1.7ffffffffffffp+0]
{1.4999999999999998})
Expect -1: Got -1 = lround(nextafter(-1.5,1) [-0x1.7ffffffffffffp+0]
{-1.4999999999999998})
Expect  1: Got  1 = lround(0.5 [0x1.0000000000000p-1] {0.5000000000000000})
Expect  2: Got  2 = lround(1.5 [0x1.8000000000000p+0] {1.5000000000000000})
Expect  3: Got  3 = lround(2.5 [0x1.4000000000000p+1] {2.5000000000000000})
Expect  4: Got  4 = lround(3.5 [0x1.c000000000000p+1] {3.5000000000000000})
Expect  0: Got  0 = lround(0.250000000000000 [0x1.fffffffffffffp-3]
{0.2500000000000000})
Expect  0: Got  0 = lround(-0.250000000000000 [-0x1.fffffffffffffp-3]
{-0.2500000000000000})
Expect  0: Got  0 = lround(0.4999999999999999 [0x1.fffffffffffffp-2]
{0.4999999999999999})
Expect  0: Got  0 = lround(-0.4999999999999999 [-0x1.fffffffffffffp-2]
{-0.4999999999999999})
Expect  1: Got  1 = lround(0.9999999999999999 [0x1.fffffffffffffp-1]
{0.9999999999999999})
Expect -1: Got -1 = lround(-0.9999999999999999 [-0x1.fffffffffffffp-1]
{-0.9999999999999999})
Expect  2: Got  2 = lround(1.9999999999999999 [0x1.fffffffffffffp+0]
{1.9999999999999998})
Expect -2: Got -2 = lround(-1.9999999999999999 [-0x1.fffffffffffffp+0]
{-1.9999999999999998})
Expect  4: Got  4 = lround(3.999999999999999 [0x1.fffffffffffffp+1]
{3.9999999999999996})
Expect -4: Got -4 = lround(-3.999999999999999 [-0x1.fffffffffffffp+1]
{-3.9999999999999996})
Expect  8: Got  8 = lround(7.999999999999999 [0x1.fffffffffffffp+2]
{7.9999999999999991})
Expect -8: Got -8 = lround(-7.999999999999999 [-0x1.fffffffffffffp+2]
{-7.9999999999999991})


-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[rsa at us dot ibm dot com]

------- Additional Comments From rsa at us dot ibm dot com  2007-11-09 19:41 -------
Buggy llround output to be used for baseline for a fix:

Expect  4503599627370496: Got   004503599627370496 = llround(2^52
[0x1.0000000000000p+52] {4503599627370496.0000000000000000})
Expect 4503599627370497 : Got 4503599627370498 = llround(2^52+1
[0x1.0000000000001p+52] {4503599627370497.0})
Expect 4503599627370498 : Got 4503599627370498 = llround(2^52+2
[0x1.0000000000002p+52] {4503599627370498.0})
Expect 4503599627370499 : Got 4503599627370500 = llround(2^52+3
[0x1.0000000000003p+52] {4503599627370499.0})
Expect 9007199254740991: Got 9007199254740992 = llround(2^52
[0x1.fffffffffffffp+52] {9007199254740991.0})
Expect 9007199254740992: Got 9007199254740992 = llround(2^52
[0x1.0000000000000p+53] {9007199254740992.0})

I have a fix for these cases that I'm trying out.

-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[rsa at us dot ibm dot com]

------- Additional Comments From rsa at us dot ibm dot com  2007-11-09 22:30 -------
The following asm stub passes all my tests cases for llround.

        lis 9,.LC0@ha
        lfd 9,.LC0@l(9)         /* Load 2^52 into fpr9.  */
        lfs 10,.LC0@l+8(9)      /* Load 0.5 into fpr10 using r9 as the base
register.  */
        fabs    2, 1            /* Get the absolute value of x.  */
        fsub    12,10,10        /* Compute 0.0 and put it in fpr12.  */
        fcmpu   6, 2, 10        /* if |x| < 0.5  */
        fcmpu   4, 2, 9         /* if |x| >= 2^52  */
        fcmpu   3, 1, 12        /* x is negative? x < 0.0  */
        blt-    6,.Lretzero     /* 0.5 > x < -0.5 so just return 0.  */
        bge-    4,.Lnobias      /* 2^52 > x < -2^52 just convert with no bias.  */
        fadd    3,2,10          /* |x|+=0.5 bias to prepare to round.  */
        bge     3,.Lconvert     /* x is positive so don't negate x.  */
        fnabs   3,3             /* -(|x|+=0.5)  */
.Lconvert:
        fctidz  4,3             /* Convert to Integer double word round toward
0.  */
        stfd    4,8(1)
        nop
        nop
        nop
        lwz     4,12(1)         /* Load return as integer.  */
        lwz     3,8(1)
         lwz 3,8(1)
.Lout:
        addi    1,1,16
        blr
.Lretzero:                      /* 0.5 > x > -0.5  */
        li      3,0             /* return 0.  */
        li      4,0
        b       .Lout
.Lnobias:
        fmr     3,2
        b       .Lconvert


-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[rsa at us dot ibm dot com]

------- Additional Comments From rsa at us dot ibm dot com  2007-11-20 21:46 -------
Created an attachment (id=2099)
 --> (http://sourceware.org/bugzilla/attachment.cgi?id=2099&action=view)
Patch to fix the erroneous lround results

This patch will fix the erroneous results in lround, lroundf, llround, and
llroundf.  It also adds some test cases to the libm-test.inc framework to
verify that they're fixed.

I want to hold back on checking this in until I get a chance to build and 'make
check' across several different systems to verify that I caught all the cases.

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
Attachment #2072 is|0                           |1
           obsolete|                            |


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[rsa at us dot ibm dot com]

------- Additional Comments From rsa at us dot ibm dot com  2007-11-20 21:50 -------
Created an attachment (id=2100)
 --> (http://sourceware.org/bugzilla/attachment.cgi?id=2100&action=view)
patch that applies -p1

This patch applies -p1

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
Attachment #2099 is|0                           |1
           obsolete|                            |


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[rsa at us dot ibm dot com]

------- Additional Comments From rsa at us dot ibm dot com  2007-11-29 16:01 -------
Patch tested on Power5, Power6 (which didn't experience the faulty results to
begin with due to extended insns for alternate rounding), PowerPC 970.  The new
make-check tests pass on all three systems.  I'm going to regenerate the patch
for CVS head.

-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[rsa at us dot ibm dot com]

------- Additional Comments From rsa at us dot ibm dot com  2007-11-29 20:59 -------
Submitted Patch to libc-alpha:
http://sources.redhat.com/ml/libc-alpha/2007-11/msg00060.html

-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[drepper at redhat dot com]

------- Additional Comments From drepper at redhat dot com  2008-04-10 05:50 -------
(In reply to comment #11)
> Created an attachment (id=2100)
 --> (http://sourceware.org/bugzilla/attachment.cgi?id=2100&action=view)
> patch that applies -p1

Doesn't work.  sysdeps/powerpc/powerpc32/power4/fpu/s_llround.S define
__llroundf which the appropriate s_llroundf.os file also defines the symbol.  If
the implementations are really identical then create an empty s_llroundf.c file
in the appropriate subdir.

And test your patches.  I see the problem on a G5.


-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |drepper at redhat dot com


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[rsa at us dot ibm dot com]

------- Additional Comments From rsa at us dot ibm dot com  2008-04-11 15:55 -------
Ulrich,

I just retried both versions of the patch (the one in the BZ and the one on the
mailing list) with both --with-cpu=power4 and --with-cpu=970 and both completed
building successfully and passed make check.

This patch didn't add or remove any files.  There should already be dummy
versions of s_llroundf.S in the necessary directories.

> find . -name s_llroundf.S
./powerpc/powerpc32/power4/fpu/s_llroundf.S
./powerpc/powerpc32/power5+/fpu/s_llroundf.S
./powerpc/powerpc32/power6/fpu/s_llroundf.S
./powerpc/powerpc64/fpu/s_llroundf.S

Did the patch apply properly?  I ask this because as an experiment I tried
removing sysdeps/powerpc/powerpc32/power4/fpu/s_llroundf.S and I _then_ got a
build warning about multiple definitions.

/home/ryanarn/build/glibc32_ppc970/math/libm_pic.a(s_llroundf.os): In function
`__llroundf':
/home/ryanarn/glibc/math/../sysdeps/powerpc/fpu/s_llroundf.c:30: multiple
definition of `__llroundf'
/home/ryanarn/build/glibc32_ppc970/math/libm_pic.a(s_llround.os):(.text+0x0):
first defined here



-- 


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.

[Bug math/4997] lround(nexafter(0.5,-1)) has incorrect value

[drepper at redhat dot com]

------- Additional Comments From drepper at redhat dot com  2008-04-11 19:33 -------
I don't know why I saw problems.  Checked in now.

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|ASSIGNED                    |RESOLVED
         Resolution|                            |FIXED


http://sourceware.org/bugzilla/show_bug.cgi?id=4997

------- You are receiving this mail because: -------
You are on the CC list for the bug, or are watching someone who is.