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

[Szabolcs Nagy <nsz@port70.net>]

* 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;
>  	}
> --