[musl] Wrong rounding in printf when precision is not set to max

[musl] Wrong rounding in printf when precision is not set to max

[John M]

I noticed printf rounding wrongly when FPU Control Word[PC] != 11.

I do not think my workaround is a fix, it just serves to show where the
reduced precision breaks the rounding.

Do you consider this a bug?
If yes, do you have an idea or comments on how a real fix would look
like?

Test case:
---
// x86_64-pc-linux-gnu-gcc float.c -o float-glibc
// x86_64-pc-linux-musl-gcc float.c -o float-musl

#include <stdio.h>

#define STREFLOP_FSTCW(cw) do { asm volatile ("fstcw %0" : "=m" (cw) : ); } while (0)
#define STREFLOP_FLDCW(cw) do { asm volatile ("fclex \n fldcw %0" : : "m" (cw)); } while (0)

void cw_set(unsigned short x87_mode) {
    STREFLOP_FLDCW(x87_mode);
}

void cw_print() {
    unsigned short x87_mode;
    STREFLOP_FSTCW(x87_mode);
    printf("FPU Control Word: 0x%04X\n", x87_mode);
}

int main() {
    cw_print();
    printf("%11.3f\n", 1200.12345);
    printf("%11.3f\n", 3.14159274);

    cw_set(0x007F);
    cw_print();
    // FIXME: musl rounds this up to 1200.124 when FPU Control Word[PC] != 11. glibc does round correctly.
    printf("%11.3f\n", 1200.12345);
    // FIXME: musl rounds this down to 3.141 when FPU Control Word[PC] != 11. glibc does round correctly.
    printf("%11.3f\n", 3.14159274);

    return 0;
}
--

My workaround:
---
diff --git a/src/stdio/vfprintf.c b/src/stdio/vfprintf.c
index 497c5e19..9b2bdb8c 100644
--- a/src/stdio/vfprintf.c
+++ b/src/stdio/vfprintf.c
@@ -13,6 +13,26 @@

 /* Some useful macros */

+#define STREFLOP_FSTCW(cw) do { __asm__ volatile ("fstcw %0" : "=m" (cw) : ); } while (0)
+#define STREFLOP_FLDCW(cw) do { __asm__ volatile ("fclex \n fldcw %0" : : "m" (cw)); } while (0)
+
+unsigned short x87_mode_old;
+void cw_push(unsigned short flags) {
+    STREFLOP_FSTCW(x87_mode_old);
+    unsigned short x87_mode = x87_mode_old | flags;
+    STREFLOP_FLDCW(x87_mode);
+}
+
+void cw_pop() {
+    STREFLOP_FLDCW(x87_mode_old);
+}
+
+void cw_print() {
+    unsigned short x87_mode;
+    STREFLOP_FSTCW(x87_mode);
+    printf("FPU Control Word: 0x%04X\n", x87_mode);
+}
+
 #define MAX(a,b) ((a)>(b) ? (a) : (b))
 #define MIN(a,b) ((a)<(b) ? (a) : (b))

@@ -318,6 +338,7 @@ static int fmt_fp(FILE *f, long double y, int w, int p, int fl, int t)
 		x = *d % i;
 		/* Are there any significant digits past j? */
 		if (x || d+1!=z) {
+			cw_push(0x0300);
 			long double round = 2/LDBL_EPSILON;
 			long double small;
 			if ((*d/i & 1) || (i==1000000000 && d>a && (d[-1]&1)))
@@ -337,6 +358,7 @@ static int fmt_fp(FILE *f, long double y, int w, int p, int fl, int t)
 				}
 				for (i=10, e=9*(r-a); *a>=i; i*=10, e++);
 			}
+			cw_pop();
 		}
 		if (z>d+1) z=d+1;
 	}
--

Re: [musl] Wrong rounding in printf when precision is not set to max

[Szabolcs Nagy]

* John M <johnm@braincalibration.de> [2024-01-01 10:11:09 +0100]:
> I noticed printf rounding wrongly when FPU Control Word[PC] != 11.
> 
> I do not think my workaround is a fix, it just serves to show where the
> reduced precision breaks the rounding.
> 
> Do you consider this a bug?

the precision control setting of x87 on i386 and x86-64 is ABI
and on linux must be double extended (64bit pc=0b11, not 53bit
pc=0b10 nor 24bit pc=0b00)

if you change the setting then the behaviour is undefined, it's
not just the libc that may misbeahave, but the compiler and
other libraries too.

> If yes, do you have an idea or comments on how a real fix would look
> like?

any x87 fp arithmetic will round the result to sinle prec float,
so the only ways to "fix" this is to use pure integer arithmetics
or change the precision setting around the internal fp computations.

i suspect glibc uses int arithmetics to compute the digits that's
why it works there, but that's just luck: glibc does not support
non-default precision setting when you call its apis either.

musl could do the same, but you would have to ensure the code
works for all supported long double formats which i guess would be
harder to do with pure int arithmetics.

> 
> Test case:
> ---
> // x86_64-pc-linux-gnu-gcc float.c -o float-glibc
> // x86_64-pc-linux-musl-gcc float.c -o float-musl
> 
> #include <stdio.h>
> 
> #define STREFLOP_FSTCW(cw) do { asm volatile ("fstcw %0" : "=m" (cw) : ); } while (0)
> #define STREFLOP_FLDCW(cw) do { asm volatile ("fclex \n fldcw %0" : : "m" (cw)); } while (0)
> 
> void cw_set(unsigned short x87_mode) {
>     STREFLOP_FLDCW(x87_mode);
> }
> 
> void cw_print() {
>     unsigned short x87_mode;
>     STREFLOP_FSTCW(x87_mode);
>     printf("FPU Control Word: 0x%04X\n", x87_mode);
> }
> 
> int main() {
>     cw_print();
>     printf("%11.3f\n", 1200.12345);
>     printf("%11.3f\n", 3.14159274);
> 
>     cw_set(0x007F);
>     cw_print();
>     // FIXME: musl rounds this up to 1200.124 when FPU Control Word[PC] != 11. glibc does round correctly.
>     printf("%11.3f\n", 1200.12345);
>     // FIXME: musl rounds this down to 3.141 when FPU Control Word[PC] != 11. glibc does round correctly.
>     printf("%11.3f\n", 3.14159274);
> 
>     return 0;
> }
> --
> 
> My workaround:
> ---
> diff --git a/src/stdio/vfprintf.c b/src/stdio/vfprintf.c
> index 497c5e19..9b2bdb8c 100644
> --- a/src/stdio/vfprintf.c
> +++ b/src/stdio/vfprintf.c
> @@ -13,6 +13,26 @@
> 
>  /* Some useful macros */
> 
> +#define STREFLOP_FSTCW(cw) do { __asm__ volatile ("fstcw %0" : "=m" (cw) : ); } while (0)
> +#define STREFLOP_FLDCW(cw) do { __asm__ volatile ("fclex \n fldcw %0" : : "m" (cw)); } while (0)
> +
> +unsigned short x87_mode_old;
> +void cw_push(unsigned short flags) {
> +    STREFLOP_FSTCW(x87_mode_old);
> +    unsigned short x87_mode = x87_mode_old | flags;
> +    STREFLOP_FLDCW(x87_mode);
> +}
> +
> +void cw_pop() {
> +    STREFLOP_FLDCW(x87_mode_old);
> +}
> +
> +void cw_print() {
> +    unsigned short x87_mode;
> +    STREFLOP_FSTCW(x87_mode);
> +    printf("FPU Control Word: 0x%04X\n", x87_mode);
> +}
> +
>  #define MAX(a,b) ((a)>(b) ? (a) : (b))
>  #define MIN(a,b) ((a)<(b) ? (a) : (b))
> 
> @@ -318,6 +338,7 @@ static int fmt_fp(FILE *f, long double y, int w, int p, int fl, int t)
>  		x = *d % i;
>  		/* Are there any significant digits past j? */
>  		if (x || d+1!=z) {
> +			cw_push(0x0300);
>  			long double round = 2/LDBL_EPSILON;
>  			long double small;
>  			if ((*d/i & 1) || (i==1000000000 && d>a && (d[-1]&1)))
> @@ -337,6 +358,7 @@ static int fmt_fp(FILE *f, long double y, int w, int p, int fl, int t)
>  				}
>  				for (i=10, e=9*(r-a); *a>=i; i*=10, e++);
>  			}
> +			cw_pop();
>  		}
>  		if (z>d+1) z=d+1;
>  	}
> --

Re: [musl] Wrong rounding in printf when precision is not set to max

[Rich Felker]

On Mon, Jan 01, 2024 at 02:43:09PM +0100, Szabolcs Nagy wrote:
> * John M <johnm@braincalibration.de> [2024-01-01 10:11:09 +0100]:
> > I noticed printf rounding wrongly when FPU Control Word[PC] != 11.
> > 
> > I do not think my workaround is a fix, it just serves to show where the
> > reduced precision breaks the rounding.
> > 
> > Do you consider this a bug?
> 
> the precision control setting of x87 on i386 and x86-64 is ABI
> and on linux must be double extended (64bit pc=0b11, not 53bit
> pc=0b10 nor 24bit pc=0b00)
> 
> if you change the setting then the behaviour is undefined, it's
> not just the libc that may misbeahave, but the compiler and
> other libraries too.

Exactly. It's UB to call any standard library code (regardless of
whether you expect it uses floating point or not) if the fpu control
word is set to a non-ABI-conforming value. If anything wants to use
alternate modes that don't match the ABI, it needs to
save/reset/restore the value around every point at which it uses
library functions.

Note that, even then, there are lots of gotchas. The x87's "single"
and "double" modes are not actually IEEE single and double, but
nonstandard types with the IEEE significand precision but excess
exponent range that's truncated via double-rounding whenever the
intermediate is spilled to memory. This will induce UB-like runaway
wrong behavior in code generated by a compiler that's not aware of
these issues, since it will effectively be able to prove 1.0==0.0. So
it is probably strongly advisable never to set these modes at all.

Rich

Re: [musl] Wrong rounding in printf when precision is not set to max

[Florian Weimer]

* Rich Felker:

> Note that, even then, there are lots of gotchas. The x87's "single"
> and "double" modes are not actually IEEE single and double, but
> nonstandard types with the IEEE significand precision but excess
> exponent range that's truncated via double-rounding whenever the
> intermediate is spilled to memory.

Doesn't setting the control word avoid double rounding?

Thanks,
Florian

Re: [musl] Wrong rounding in printf when precision is not set to max

[Rich Felker]

On Tue, Jan 02, 2024 at 11:16:15AM +0100, Florian Weimer wrote:
> * Rich Felker:
> 
> > Note that, even then, there are lots of gotchas. The x87's "single"
> > and "double" modes are not actually IEEE single and double, but
> > nonstandard types with the IEEE significand precision but excess
> > exponent range that's truncated via double-rounding whenever the
> > intermediate is spilled to memory.
> 
> Doesn't setting the control word avoid double rounding?

No. The first rounding is of the exact mathematical result to
x87-single (23-bit significand, 16-bit exponent) or x87-double (52-bit
significand, 16-bit exponent), then there's a second rounding on store
to IEEE single or IEEE double. The result differs on overflow or
underflow. For example:

- a+b-c (including special case a+b-b) can be finite in x87-double
  (even after final store), but infinite in IEEE double.

- a+b-c (including special case a+b-b) can be exact in x87-double, but
  incur rounding in IEEE double due to a+b becoming denormal.

If you can be certain no intermediates will ever overflow or underflow
(and no denormals appear in inputs), then I believe x87 reduced
precision modes can faithfully emulate IEEE arithmetic. But this is a
nontrivial condition to satisfy in general.

Rich