Follow @Openwall on Twitter for new release announcements and other news
[<prev] [next>] [<thread-prev] [thread-next>] [day] [month] [year] [list]
Message-ID: <20250701161956.GC1827@brightrain.aerifal.cx>
Date: Tue, 1 Jul 2025 12:19:57 -0400
From: Rich Felker <dalias@...c.org>
To: Michal Terepeta <michalt@...gle.com>
Cc: musl@...ts.openwall.com, t5y-external <t5y-external@...gle.com>
Subject: Re: Wrong formatting of doubles?

On Tue, Jul 01, 2025 at 10:22:25AM -0400, Rich Felker wrote:
> On Tue, Jul 01, 2025 at 07:55:20AM +0200, Michal Terepeta wrote:
> > Hi,
> > 
> > We're working on a RISC-V system and using musl as our libc. A recent change
> > in musl seems to have caused wrong formatting (printf) of double values:
> > https://git.musl-libc.org/cgit/musl/commit/src/stdio/vfprintf.c?id=572a2e2eb91f00f2f25d301cfb50f435e7ae16b3
> > 
> > Using a simple binary to print max double [1] with the current code I get:
> > 
> > ```
> > The maximum value of a double (printf %e): 1.348676e+308
> > The maximum value of a double (printf %E): 1.348676E+308
> > The maximum value of a double (printf %g): 3.13486e+308
> > ```
> > 
> > With the patch [2] that reverts most of the above change, I get the expected
> > output:
> > 
> > ```
> > The maximum value of a double (printf %e): 1.797693e+308
> > The maximum value of a double (printf %E): 1.797693E+308
> > The maximum value of a double (printf %g): 1.79769e+308
> > ```
> > 
> > It'd be great if someone could take a look if they can also repro and perhaps
> > revert the change?
> > 
> > Many thanks!
> > 
> > Michal
> > 
> > 
> > 
> > [1] Repro program:
> > 
> > ```
> > #include <float.h>
> > #include <stdio.h>
> > 
> > int main() {
> >   printf("The maximum value of a double (printf %%e): %e\n", DBL_MAX);
> >   printf("The maximum value of a double (printf %%E): %E\n", DBL_MAX);
> >   printf("The maximum value of a double (printf %%g): %g\n", DBL_MAX);
> >   return 0;
> > }
> > 
> > ```
> > 
> > 
> > [2] Patch to test:
> > 
> > ```
> > diff --git a/src/stdio/vfprintf.c b/src/stdio/vfprintf.c
> > index 010041ca6c..01004158bf 100644
> > --- a/src/stdio/vfprintf.c
> > +++ b/src/stdio/vfprintf.c
> > @@ -180,12 +180,8 @@
> > 
> >  static int fmt_fp(FILE *f, long double y, int w, int p, int fl, int t, int ps)
> >  {
> > -       int bufsize = (ps==BIGLPRE)
> > -               ? (LDBL_MANT_DIG+28)/29 + 1 +          // mantissa expansion
> > -                 (LDBL_MAX_EXP+LDBL_MANT_DIG+28+8)/9  // exponent expansion
> > -               : (DBL_MANT_DIG+28)/29 + 1 +
> > -                 (DBL_MAX_EXP+DBL_MANT_DIG+28+8)/9;
> > -       uint32_t big[bufsize];
> > +       uint32_t big[(LDBL_MANT_DIG+28)/29 + 1          // mantissa expansion
> > +               + (LDBL_MAX_EXP+LDBL_MANT_DIG+28+8)/9]; // exponent expansion
> >         uint32_t *a, *d, *r, *z;
> >         int e2=0, e, i, j, l;
> >         char buf[9+LDBL_MANT_DIG/4], *s;
> > ```
> 
> Could you try the attached patch?
> 

> diff --git a/src/stdio/vfprintf.c b/src/stdio/vfprintf.c
> index 76733997..2d2f4f3e 100644
> --- a/src/stdio/vfprintf.c
> +++ b/src/stdio/vfprintf.c
> @@ -180,11 +180,10 @@ typedef char compiler_defines_long_double_incorrectly[9-(int)sizeof(long double)
>  
>  static int fmt_fp(FILE *f, long double y, int w, int p, int fl, int t, int ps)
>  {
> -	int bufsize = (ps==BIGLPRE)
> -		? (LDBL_MANT_DIG+28)/29 + 1 +          // mantissa expansion
> -		  (LDBL_MAX_EXP+LDBL_MANT_DIG+28+8)/9  // exponent expansion
> -		: (DBL_MANT_DIG+28)/29 + 1 +
> -		  (DBL_MAX_EXP+DBL_MANT_DIG+28+8)/9;
> +	int max_mant_dig = (ps==BIGLPRE) ? LDBL_MANT_DIG : DBL_MANT_DIG;
> +	int max_exp = (ps==BIGLPRE) ? LDBL_MAX_EXP : DBL_MAX_EXP;
> +	int bufsize = (max_mant_dig+28)/29 + 1          // mantissa expansion
> +	            + (max_exp+max_mant_dig+28+8)/9;    // exponent expansion
>  	uint32_t big[bufsize];
>  	uint32_t *a, *d, *r, *z;
>  	int e2=0, e, i, j, l;
> @@ -266,7 +265,7 @@ static int fmt_fp(FILE *f, long double y, int w, int p, int fl, int t, int ps)
>  	if (y) y *= 0x1p28, e2-=28;
>  
>  	if (e2<0) a=r=z=big;
> -	else a=r=z=big+sizeof(big)/sizeof(*big) - LDBL_MANT_DIG - 1;
> +	else a=r=z=big+sizeof(big)/sizeof(*big) - max_mant_dig - 1;
>  
>  	do {
>  		*z = y;

To clarify, the root cause here is that the subtraction later down the
function to make room for the mantissa, fixed in the last quoted hunk
above, was still using a hardcoded LDBL_MANT_DIG rather than adapting
to the size of the type being formatted.

I was able to reproduce the buf on aarch64, albeit with all-zeros
output instead of the prefixed '3'. Hitting it seems to depend on
LDBL_MANT_DIG-DBL_MANT_DIG being high enough, which only happens on
archs with quad ld (not ld80).

For what it's worth, the math does not look right even with the buffer
sized for LDBL_MANT_DIG, and it's probably only by luck that it
worked. We only reserved one b1b (base-1000000000) slot per 29 bits of
mantissa, but adjusted the end pointer to take away a whole b1b slot
for each bit of the mantissa, potentially leaving too little room for
the exponent expansion. In practice this worked because positive
exponent expansions use a lot less (roughly 1/3) the estimated space;
it's negative exponent expansions that use the full extimate. (This is
because each halving adds an extra digit to the end, but it takes
log2(10) doublings to add a digit to the front.)

I think the correct code would be something like:

+	else a=r=z=big+sizeof(big)/sizeof(*big) - (max_mant_dig+28)/29 - 1;

And maybe this math should be done once at the top and reused rather
than repeated (i.e. use max_mant_slots instead of max_mant_dig). This
would also eliminate the need for the comments by documenting the
intermediate size calculations in variable names.


So my understanding of the impact is this:

- The newly exposed bug from using a smaller buffer for double is not
  present in any release, and only affected some archs (those with
  long double as quad).

- The longstanding conceptual bug is not an actual bug, only because
  luck of the particular scales of the mantissa and exponent size
  parameters for long double avoided going over the excess space we
  have in the positive-exponent case.

- Anyone who's been running latest-from-git does need to apply a fix.



Powered by blists - more mailing lists

Confused about mailing lists and their use? Read about mailing lists on Wikipedia and check out these guidelines on proper formatting of your messages.