[Bug libc/5044] New: printf doesn't take the rounding mode into account

[Bug libc/5044] New: printf doesn't take the rounding mode into account

[vincent+libc at vinc17 dot org]

When printing out the value of a double, e.g. with %f, the rounding mode isn't
taken into account. For instance, the program below:

#include <stdio.h>
#include <fenv.h>

static volatile double x = 1.0 / 3.0;

void out (const char *s, int r)
{
  if (fesetround (r))
    fprintf (stderr, "fesetround error\n");
  else
    printf ("1/3 rounded %s: %1.6f\n", s, x);
}

int main (void)
{
  out ("downward", FE_DOWNWARD);
  out ("  upward", FE_UPWARD);
  return 0;
}

gives:

1/3 rounded downward: 0.333333
1/3 rounded   upward: 0.333333

Since (1.0/3.0) isn't exactly representable in decimal, the rounded decimal
values should be different (more precisely, they should be 0.333333 and 0.333334).

The ISO C standard says in 7.19.6.1#13:

"For e, E, f, F, g, and G conversions, if the number of significant decimal
digits is at most DECIMAL_DIG, then the result should be correctly rounded.243)
If the number of significant decimal digits is more than DECIMAL_DIG but the
source value is exactly representable with DECIMAL_DIG digits, then the result
should be an exact representation with trailing zeros. Otherwise, the source
value is bounded by two adjacent decimal strings L < U, both having DECIMAL_DIG
significant digits; the value of the resultant decimal string D should satisfy L
 <= D <= U, with the extra stipulation that the error should have a correct sign
for the current rounding direction."

-- 
           Summary: printf doesn't take the rounding mode into account
           Product: glibc
           Version: unspecified
            Status: NEW
          Severity: normal
          Priority: P2
         Component: libc
        AssignedTo: drepper at redhat dot com
        ReportedBy: vincent+libc at vinc17 dot org
                CC: glibc-bugs at sources dot redhat dot com
  GCC host triplet: i686-pc-linux-gnu, powerpc-unknown-linux-gnu


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[eberlein at us dot ibm dot com]

------- Additional Comments From eberlein at us dot ibm dot com  2008-04-15 21:52 -------
Checking the rounding mode from printf would be done using fegetround, which
would create a dependency on libm from libc.  Comments?

-- 


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[vincent+libc at vinc17 dot org]

------- Additional Comments From vincent+libc at vinc17 dot org  2008-04-15 22:53 -------
Checking the rounding mode can also be done with two well-chosen additions, for
instance. So, no dependency on libm is needed.

-- 


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[rsa at us dot ibm dot com]

------- Additional Comments From rsa at us dot ibm dot com  2008-04-16 14:20 -------
You can skirt the libm dependency by using the internal glibc macros _FPU_GETCW
and _FPU_SETCW.  On architectures which have hardware rounding mode, like x86_64
and powerpc64 it'll access the rounding mode directly in the registers (cw
register and FPSCR register respectively).

On architectures without a hardware rounding mode it'll use the rounding mode
that is already stored in the glibc thread-local-storage.

-- 


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[rsa at us dot ibm dot com]

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |eberlein at us dot ibm dot
                   |                            |com


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[rsa at us dot ibm dot com]

-- 
           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |rsa at us dot ibm dot com


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[eberlein at us dot ibm dot com]

------- Additional Comments From eberlein at us dot ibm dot com  2008-04-17 20:28 -------
Created an attachment (id=2706)
 --> (http://sourceware.org/bugzilla/attachment.cgi?id=2706&action=view)
Quick patch to allow printf_fp to get rounding mode using _FPU_GETCW


-- 


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[eberlein at us dot ibm dot com]

------- Additional Comments From eberlein at us dot ibm dot com  2008-04-17 20:33 -------
I've attached a patch based on Ryan's suggestion to use _FPU_GETCW.  It fixes
the example given in the description.  It may have some problems with the
existing ties-to-even code.  I've only tested it on POWER.  Please check that
the patch works on other architectures, and I'll work on the ties-to-even bit.

-- 


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[vincent+libc at vinc17 dot org]

------- Additional Comments From vincent+libc at vinc17 dot org  2008-04-17 23:20 -------
I haven't looked at the code itself, but the patch looks suspicious to me.
Shouldn't L'0' be replaced by L'0' - 1 so that in directed rounding modes, the
test "digit > rounddigit" is either always true or always false.

-- 


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[eberlein at us dot ibm dot com]

------- Additional Comments From eberlein at us dot ibm dot com  2008-04-17 23:42 -------
(In reply to comment #6)
> Shouldn't L'0' be replaced by L'0' - 1 so that in directed rounding modes, the
> test "digit > rounddigit" is either always true or always false.

No. You never round when the digit is zero.

-- 


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[vincent+libc at vinc17 dot org]

------- Additional Comments From vincent+libc at vinc17 dot org  2008-04-18 07:34 -------
(In reply to comment #7)
> No. You never round when the digit is zero.

Hmm, I was assuming that the exact values were taken apart... There seems to be
something missing in the code then. For instance, how does
  printf("%1.2f\n", 0.5001);
round the value upward (0.51) in the FE_UPWARD rounding mode? I think you need
some additional flag to say whether the converted value is exact or not, i.e.
whether the remaining digits are all zeros.

-- 


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[eberlein at us dot ibm dot com]

------- Additional Comments From eberlein at us dot ibm dot com  2008-04-18 20:57 -------
Created an attachment (id=2713)
 --> (http://sourceware.org/bugzilla/attachment.cgi?id=2713&action=view)
Revised patch to allow printf_fp to get rounding mode using _FPU_GETCW

Vincent, I see what you mean about taking into account whether the rounded
digit is followed by zeroes.  Rounding upward should never result in a value
that is less than the original value, which wasn't caught when the rounded
digit is zero followed by nonzero digits.  I'm attaching a patch that modifies
the existing ties-to-even code to work for FE_UPWARD and FE_DOWNWARD that
checks for trailing zeroes when the rounded digit is zero.

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


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[vincent+libc at vinc17 dot org]

------- Additional Comments From vincent+libc at vinc17 dot org  2008-04-18 23:58 -------
(In reply to comment #10)
> For example 0.5001 is really 0x1.000d1b71758e20p-1

In my example, the exact value doesn't matter: the rounding to two digits will
be the same for any good approximation to 0.5001 (denoted ~0.5001 below).

> "Otherwise, the source value is bounded by two adjacent decimal strings L < U,
> both having DECIMAL_DIG significant digits; the value of the resultant decimal
> string D should satisfy L <= D <= U, with the extra stipulation that the error
> should have a correct sign for the current rounding direction."
> 
> This implies that either L or U (.50 or .51) are allowed!

No: if .50 is output, the error does not have a correct sign in rounding upward,
since 0.50 < ~0.5001. So, only .51 is allowed here. Now, this part of the ISO C
standard does not apply anyway since 2 digits are required and 2 <= DECIMAL_DIG.
You need to look at the beginning of the paragraph:

"For e, E, f, F, g, and G conversions, if the number of significant decimal
digits is at most DECIMAL_DIG, then the result should be correctly rounded."

and the correctly-rounded value is .51 here.

-- 


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[sjmunroe at us dot ibm dot com]

------- Additional Comments From sjmunroe at us dot ibm dot com  2008-04-18 22:03 -------
the core problem is that you are trying to apply decimal rounding rules to
IEEE754 Binary float.

For example 0.5001 is really 0x1.000d1b71758e20p-1 (format %.14a which shows the
mantisa in hex). Note that there is no exact binary representation for rounding
value 0.005 implied by this request.

To get what Vincent says he wants would require converting the IEEE754 double to
decimal and perform the rounding in decimal (via quantize). It is is not clear
that The ISO C standard says in 7.19.6.1#13 requires this:

"Otherwise, the source value is bounded by two adjacent decimal strings L < U,
both having DECIMAL_DIG significant digits; the value of the resultant decimal
string D should satisfy L <= D <= U, with the extra stipulation that the error
should have a correct sign for the current rounding direction."

This implies that either L or U (.50 or .51) are allowed!

-- 


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

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

[Bug libc/5044] printf doesn't take the rounding mode into account

[khalil dot ghorbal at cea dot fr]

------- Additional Comments From khalil dot ghorbal at cea dot fr  2008-06-05 17:52 -------
Hi all,

I have noticed the same problem with libc-2.3.6 (I'm using a precompiled libc6
debian(etch) package).

In the example submitted by Vincent, maybe we should calculate x after choosing
the rounding mode: since x
is evaluated first, the default rounding mode (toward +oo) is used. Altering the
rounding mode later won't
change the internal binary representation for x (the type qualifier volatile
here is without effect).

But even doing this don't change the output of printf !

In fact, to get a correct (with respect to IEEE 754) result, I used an
intermediate variable, say "b", set
to 3.0, we get (see bottom for modified source code):
(using %a)
1/3 rounded downward: 0x1.5555555555555p-2
1/3 rounded   upward: 0x1.5555555555556p-2
(using %.17f)
1/3 rounded downward: 0.33333333333333331
1/3 rounded   upward: 0.33333333333333337

A height precision (at least .17) is needed to see the difference using %f
instead of %a (and this is not a normal behaviour, as one should see directly
the difference, since we print two different numbers !).

To summarize :
- fesetround don't influence the internal representation of constants such as
1./3. (or 1./10.)
- using %f seems to have some troubles ...


Hope this helps.

Regards,

<<<<<< source modified >>>>>>>>
#include <stdio.h>
#include <fenv.h>

static volatile double x;

void out (const char *s, int r)
{
        double b = 3.0;
        if (fesetround (r))
                fprintf (stderr, "fesetround error\n");
        else
                x = 1.0/b;
        printf ("1/3 rounded %s: %a\n", s, x);
}

int main (void)
{
        out ("downward", FE_DOWNWARD);
        out ("  upward", FE_UPWARD);
        return 0;
}

-- 


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

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