Bug 5044 - printf doesn't take the rounding mode into account
Summary: printf doesn't take the rounding mode into account
Status: RESOLVED FIXED
Alias: None
Product: glibc
Classification: Unclassified
Component: stdio (show other bugs)
Version: unspecified
: P2 normal
Target Milestone: ---
Assignee: Ulrich Drepper
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2007-09-19 13:44 UTC by Vincent Lefèvre
Modified: 2014-07-04 07:29 UTC (History)
3 users (show)

See Also:
Host: i686-pc-linux-gnu, powerpc-unknown-linux-gnu
Target:
Build:
Last reconfirmed:
Project(s) to access:
ssh public key:
fweimer: security-


Attachments
Quick patch to allow printf_fp to get rounding mode using _FPU_GETCW (648 bytes, patch)
2008-04-17 20:28 UTC, Pete Eberlein
Details | Diff
Revised patch to allow printf_fp to get rounding mode using _FPU_GETCW (745 bytes, patch)
2008-04-18 20:57 UTC, Pete Eberlein
Details | Diff

Note You need to log in before you can comment on or make changes to this bug.
Description Vincent Lefèvre 2007-09-19 13:44:19 UTC
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."
Comment 1 Pete Eberlein 2008-04-15 21:52:08 UTC
Checking the rounding mode from printf would be done using fegetround, which
would create a dependency on libm from libc.  Comments?
Comment 2 Vincent Lefèvre 2008-04-15 22:53:35 UTC
Checking the rounding mode can also be done with two well-chosen additions, for
instance. So, no dependency on libm is needed.
Comment 3 Ryan S. Arnold 2008-04-16 14:20:24 UTC
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.
Comment 4 Pete Eberlein 2008-04-17 20:28:05 UTC
Created attachment 2706 [details]
Quick patch to allow printf_fp to get rounding mode using _FPU_GETCW
Comment 5 Pete Eberlein 2008-04-17 20:33:12 UTC
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.
Comment 6 Vincent Lefèvre 2008-04-17 23:20:59 UTC
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.
Comment 7 Pete Eberlein 2008-04-17 23:42:16 UTC
(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.
Comment 8 Vincent Lefèvre 2008-04-18 07:34:14 UTC
(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.
Comment 9 Pete Eberlein 2008-04-18 20:57:30 UTC
Created attachment 2713 [details]
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.
Comment 10 Steven Munroe 2008-04-18 22:03:52 UTC
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!
Comment 11 Vincent Lefèvre 2008-04-18 23:58:09 UTC
(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.
Comment 12 Khalil Ghorbal 2008-06-05 17:52:32 UTC
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;
}
Comment 13 Joseph Myers 2012-09-14 20:24:24 UTC
The case of decimal output (%e/%E/%f/%F/%g/%G) is fixed for 2.17 by:

commit 784761bee3828e4e21cbe15340c7a99e3fd1788b
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Fri Sep 14 20:18:49 2012 +0000

    Make printf respect the rounding mode for decimal output (bug 5044).

That is, decimal output should now be correctly rounding for all rounding modes, where previously it was for round-to-nearest only.

Hex output (%a/%A) has the same issue and still needs to be fixed.
Comment 14 Joseph Myers 2012-09-24 15:40:18 UTC
The remaining case, that of hex output (%a, %A), is fixed for 2.17 by:

commit a9f8e53a5b14ba481999ded036b025554ea06362
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Mon Sep 24 15:38:21 2012 +0000

    Make printf respect the rounding mode for hex output (bug 5044).
Comment 15 Jackie Rosen 2014-02-16 18:25:34 UTC Comment hidden (spam)