Wrong formatting of doubles?

Wrong formatting of doubles?

[Michal Terepeta]

[-- Attachment #1: Type: text/plain, Size: 2075 bytes --]

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;
```

[-- Attachment #2: S/MIME Cryptographic Signature --]
[-- Type: application/pkcs7-signature, Size: 5279 bytes --]

Re: [musl] Wrong formatting of doubles?

[Rich Felker]

[-- Attachment #1: Type: text/plain, Size: 2318 bytes --]

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?


[-- Attachment #2: printf_dbl_bug.diff --]
[-- Type: text/plain, Size: 1207 bytes --]

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;

Re: Wrong formatting of doubles?

[Rich Felker]

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.

Re: Wrong formatting of doubles?

[Rich Felker]

On Tue, Jul 01, 2025 at 12:19:57PM -0400, Rich Felker wrote:
> 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;

An additional -1 (+1 to the # of slots for mantissa expansion, as in
the commented expression) is needed here because the do/while loop
emits a final zero slot for the mantissa. I'm not sure this is
actually needed later, but as long as it's there it needs to be
accounted for.

And indeed, with some debug instrumentation, empirically
z==buf+bufsize after the initial mantissa expansion loop finishes.

> 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.

As a result of the above findings, I'm strongly in favor of doing it
this way. Repeating a confusing expression like this is a formula for
disaster.

I'd also like to think about ways we could avoid introducing bugs like
this in the future. In this case, just a testcase for printing DBL_MAX
in multiple formats would likely have caught it. Simple trapping-only
UBsan might have caught it, but I'm not sure if that would follow all
the positional arithmetic in fmt_fp to detect going outside the object
bounds. This may be worth investigating.

Big thanks to Michal Terepeta for catching this prior to it slipping
into a release!

Rich

Re: Wrong formatting of doubles?

[Jeffrey Walton]

On Tue, Jul 1, 2025 at 12:38 PM Rich Felker <dalias@libc.org> 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
> > > > [...]
> [...]
> I'd also like to think about ways we could avoid introducing bugs like
> this in the future. In this case, just a testcase for printing DBL_MAX
> in multiple formats would likely have caught it. Simple trapping-only
> UBsan might have caught it, but I'm not sure if that would follow all
> the positional arithmetic in fmt_fp to detect going outside the object
> bounds. This may be worth investigating.

Yes, test cases under a tool like a memory checker or UB checker are
the way to go. In this case, the sanitizers or Valgrind likely would
have alerted to the problematic code.

Another trick I like to use in debug builds is the use of asserts to
verify all program state. In this case, you might be able to check if
a pointer is within a particular range. In release builds, the checks
can be omitted since they were already verified using a test case in
debug builds, so the code will not take a performance hit.

When I perform root cause analysis on why a bug hit production, it
just about always (always?) includes a missing test case.

Jeff

Re: [musl] Wrong formatting of doubles?

[Markus Wichmann]

Am Tue, Jul 01, 2025 at 12:37:58PM -0400 schrieb Rich Felker:
> On Tue, Jul 01, 2025 at 12:19:57PM -0400, Rich Felker wrote:
> > On Tue, Jul 01, 2025 at 10:22:25AM -0400, Rich Felker wrote:
> > >  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;

I am wondering how this is causing wrong digits instead of a crash. The
only thing the above changes could prevent is an out-of-bounds access.

Well, OK, let's do some math:
DBL_MANT_DIG = 53
LDBL_MANT_DIG = 113
DBL_MAX_EXP = 1024
bufsize = 126

That only leaves 13 blocks for the exponent expansion, which is enough
for 117 decimal digits in front of the radix before a buffer overflow
occurs. Regrettably, DBL_MAX has 308 digits in front of the radix. But
that should mean that some return address gets overwritten and a crash
occurs on return, right?

With the corrected calculation, more than double the required digits can
be saved.

> > I think the correct code would be something like:
> > 
> > +	else a=r=z=big+sizeof(big)/sizeof(*big) - (max_mant_dig+28)/29 - 1;
> 
> An additional -1 (+1 to the # of slots for mantissa expansion, as in
> the commented expression) is needed here because the do/while loop
> emits a final zero slot for the mantissa. I'm not sure this is
> actually needed later, but as long as it's there it needs to be
> accounted for.
> 
> And indeed, with some debug instrumentation, empirically
> z==buf+bufsize after the initial mantissa expansion loop finishes.
> 

I am also wondering if there is any need for the calculation to be exact
in this case. The buffer size is as big as it is for negative exponents,
and for positive ones it really doesn't matter, since the decimal
expansion can never be big enough to fill the whole buffer. So the only
thing this changes is where exactly the data is inside the buffer, which
doesn't seem like a worthwhile change to me.

Ciao,
Markus

Re: Wrong formatting of doubles?

[Rich Felker]

On Tue, Jul 01, 2025 at 07:22:20PM +0200, Markus Wichmann wrote:
> Am Tue, Jul 01, 2025 at 12:37:58PM -0400 schrieb Rich Felker:
> > On Tue, Jul 01, 2025 at 12:19:57PM -0400, Rich Felker wrote:
> > > On Tue, Jul 01, 2025 at 10:22:25AM -0400, Rich Felker wrote:
> > > >  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;
> 
> I am wondering how this is causing wrong digits instead of a crash. The
> only thing the above changes could prevent is an out-of-bounds access.

Presumably it's overflowing into other locals on the stack and
clobbering them in a way that affects the subsequent flow of fmt_fp.

> Well, OK, let's do some math:
> DBL_MANT_DIG = 53
> LDBL_MANT_DIG = 113
> DBL_MAX_EXP = 1024
> bufsize = 126
> 
> That only leaves 13 blocks for the exponent expansion, which is enough
> for 117 decimal digits in front of the radix before a buffer overflow
> occurs. Regrettably, DBL_MAX has 308 digits in front of the radix. But
> that should mean that some return address gets overwritten and a crash
> occurs on return, right?

Nope. The overflow is downward (lower addresses) but the return
address and most of the static contents of the stack frame should be
upward relative the the VLA.

I'm not actually sure how things end up being below the VLA in a way
that causes the erroneous output. FWIW on aarch64 it was all-zeros.
But in any case it's clearly UB, and the mechanism isn't that
important unless you're trying to write an exploit.

> With the corrected calculation, more than double the required digits can
> be saved.
> 
> > > I think the correct code would be something like:
> > > 
> > > +	else a=r=z=big+sizeof(big)/sizeof(*big) - (max_mant_dig+28)/29 - 1;
> > 
> > An additional -1 (+1 to the # of slots for mantissa expansion, as in
> > the commented expression) is needed here because the do/while loop
> > emits a final zero slot for the mantissa. I'm not sure this is
> > actually needed later, but as long as it's there it needs to be
> > accounted for.
> > 
> > And indeed, with some debug instrumentation, empirically
> > z==buf+bufsize after the initial mantissa expansion loop finishes.
> 
> I am also wondering if there is any need for the calculation to be exact
> in this case. The buffer size is as big as it is for negative exponents,
> and for positive ones it really doesn't matter, since the decimal
> expansion can never be big enough to fill the whole buffer. So the only
> thing this changes is where exactly the data is inside the buffer, which
> doesn't seem like a worthwhile change to me.

It doesn't need to be exact, but if it's not, we're relying on roughly
28/29 of max_mant_dig being smaller than 2/27 of max_exp. In practice
this is true for the paramters of real-world floating point formats,
but it's not reasonable as an undocumented assumption to rely on for
safety. Having the values be correctly and exactly computed makes them
self-documenting.

Rich

Re: Wrong formatting of doubles?

[Rich Felker]

[-- Attachment #1: Type: text/plain, Size: 7173 bytes --]

On Tue, Jul 01, 2025 at 12:37:58PM -0400, Rich Felker wrote:
> On Tue, Jul 01, 2025 at 12:19:57PM -0400, Rich Felker wrote:
> > 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;
> 
> An additional -1 (+1 to the # of slots for mantissa expansion, as in
> the commented expression) is needed here because the do/while loop
> emits a final zero slot for the mantissa. I'm not sure this is
> actually needed later, but as long as it's there it needs to be
> accounted for.
> 
> And indeed, with some debug instrumentation, empirically
> z==buf+bufsize after the initial mantissa expansion loop finishes.
> 
> > 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.
> 
> As a result of the above findings, I'm strongly in favor of doing it
> this way. Repeating a confusing expression like this is a formula for
> disaster.
> 
> I'd also like to think about ways we could avoid introducing bugs like
> this in the future. In this case, just a testcase for printing DBL_MAX
> in multiple formats would likely have caught it. Simple trapping-only
> UBsan might have caught it, but I'm not sure if that would follow all
> the positional arithmetic in fmt_fp to detect going outside the object
> bounds. This may be worth investigating.
> 
> Big thanks to Michal Terepeta for catching this prior to it slipping
> into a release!

Proposed commit attached. Comments welcome. I'll plan to push tomorrow
if there are no objections.


[-- Attachment #2: 0001-printf-fix-regression-in-large-double-formatting-on-.patch --]
[-- Type: text/plain, Size: 2948 bytes --]

From f96e47a26102d537c29435f0abf9ec94676a030e Mon Sep 17 00:00:00 2001
From: Rich Felker <dalias@aerifal.cx>
Date: Tue, 1 Jul 2025 21:30:18 -0400
Subject: [PATCH] printf: fix regression in large double formatting on ld128
 archs

commit 572a2e2eb91f00f2f25d301cfb50f435e7ae16b3 adjusted the buffer
for decimal conversion to be a VLA that only uses the full size needed
for long double when the argument type was long double. however, it
failed to update a later expression for the positioning within the
buffer, which still used a fixed offset of LDBL_MANT_DIG. this caused
doubles with a large positive exponent to overflow below the start of
the array, producing wrong output and potentially runaway wrong
execution.

this bug has not been present in any release, and has not been
analyzed in depth for security considerations.

it turns out the original buffer offset expression involving
LDBL_MANT_DIG was incorrect as well, and only worked because the space
reserved for expanding the exponent is roughly 3 times the size it
needs to be when the exponent is positive, leaving plenty of extra
space to compensate for the error. the actual offset should be in
base-1000000000 slot units, not bits, and numerically equal to the
number of slots that were previously allocated for mantissa expansion.

in order to ensure consistency and make the code more comprehensible,
commented subexpressions are replaced by intermediate named variables,
and the newly introduced max_mant_slots is used for both the
allocation and the buffer offset adjustment. the included +1 term
accounts for a trailing zero slot that's always emitted.
---
 src/stdio/vfprintf.c | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/stdio/vfprintf.c b/src/stdio/vfprintf.c
index 76733997..a68edabb 100644
--- a/src/stdio/vfprintf.c
+++ b/src/stdio/vfprintf.c
@@ -180,11 +180,11 @@ 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 max_mant_slots = (max_mant_dig+28)/29 + 1;
+	int max_exp_slots = (max_exp+max_mant_dig+28+8)/9;
+	int bufsize = max_mant_slots + max_exp_slots;
 	uint32_t big[bufsize];
 	uint32_t *a, *d, *r, *z;
 	int e2=0, e, i, j, l;
@@ -266,7 +266,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_slots - 1;
 
 	do {
 		*z = y;
-- 
2.21.0

Re: [musl] Wrong formatting of doubles?

[Michal Terepeta]

[-- Attachment #1: Type: text/plain, Size: 328 bytes --]

On Wed, Jul 2, 2025 at 3:46 AM Rich Felker <dalias@libc.org> wrote:
> [...]
>
> Proposed commit attached. Comments welcome. I'll plan to push tomorrow
> if there are no objections.
>

I just tested the patch on our end (64-bit RISC-V) and everything looks good!
Thank you for a quick response and fix! :D

- Michal

[-- Attachment #2: S/MIME Cryptographic Signature --]
[-- Type: application/pkcs7-signature, Size: 5279 bytes --]

Re: [musl] Wrong formatting of doubles?

[Rich Felker]

On Tue, Jul 01, 2025 at 09:46:10PM -0400, Rich Felker wrote:
> On Tue, Jul 01, 2025 at 12:37:58PM -0400, Rich Felker wrote:
> > On Tue, Jul 01, 2025 at 12:19:57PM -0400, Rich Felker wrote:
> > > 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;
> > 
> > An additional -1 (+1 to the # of slots for mantissa expansion, as in
> > the commented expression) is needed here because the do/while loop
> > emits a final zero slot for the mantissa. I'm not sure this is
> > actually needed later, but as long as it's there it needs to be
> > accounted for.
> > 
> > And indeed, with some debug instrumentation, empirically
> > z==buf+bufsize after the initial mantissa expansion loop finishes.
> > 
> > > 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.
> > 
> > As a result of the above findings, I'm strongly in favor of doing it
> > this way. Repeating a confusing expression like this is a formula for
> > disaster.
> > 
> > I'd also like to think about ways we could avoid introducing bugs like
> > this in the future. In this case, just a testcase for printing DBL_MAX
> > in multiple formats would likely have caught it. Simple trapping-only
> > UBsan might have caught it, but I'm not sure if that would follow all
> > the positional arithmetic in fmt_fp to detect going outside the object
> > bounds. This may be worth investigating.
> > 
> > Big thanks to Michal Terepeta for catching this prior to it slipping
> > into a release!
> 
> Proposed commit attached. Comments welcome. I'll plan to push tomorrow
> if there are no objections.
> 

> >From f96e47a26102d537c29435f0abf9ec94676a030e Mon Sep 17 00:00:00 2001
> From: Rich Felker <dalias@aerifal.cx>
> Date: Tue, 1 Jul 2025 21:30:18 -0400
> Subject: [PATCH] printf: fix regression in large double formatting on ld128
>  archs
> 
> commit 572a2e2eb91f00f2f25d301cfb50f435e7ae16b3 adjusted the buffer
> for decimal conversion to be a VLA that only uses the full size needed
> for long double when the argument type was long double. however, it
> failed to update a later expression for the positioning within the
> buffer, which still used a fixed offset of LDBL_MANT_DIG. this caused
> doubles with a large positive exponent to overflow below the start of
> the array, producing wrong output and potentially runaway wrong
> execution.
> 
> this bug has not been present in any release, and has not been
> analyzed in depth for security considerations.
> 
> it turns out the original buffer offset expression involving
> LDBL_MANT_DIG was incorrect as well, and only worked because the space
> reserved for expanding the exponent is roughly 3 times the size it
> needs to be when the exponent is positive, leaving plenty of extra
> space to compensate for the error. the actual offset should be in
> base-1000000000 slot units, not bits, and numerically equal to the
> number of slots that were previously allocated for mantissa expansion.
> 
> in order to ensure consistency and make the code more comprehensible,
> commented subexpressions are replaced by intermediate named variables,
> and the newly introduced max_mant_slots is used for both the
> allocation and the buffer offset adjustment. the included +1 term
> accounts for a trailing zero slot that's always emitted.
> ---
>  src/stdio/vfprintf.c | 12 ++++++------
>  1 file changed, 6 insertions(+), 6 deletions(-)
> 
> diff --git a/src/stdio/vfprintf.c b/src/stdio/vfprintf.c
> index 76733997..a68edabb 100644
> --- a/src/stdio/vfprintf.c
> +++ b/src/stdio/vfprintf.c
> @@ -180,11 +180,11 @@ 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 max_mant_slots = (max_mant_dig+28)/29 + 1;
> +	int max_exp_slots = (max_exp+max_mant_dig+28+8)/9;
> +	int bufsize = max_mant_slots + max_exp_slots;
>  	uint32_t big[bufsize];
>  	uint32_t *a, *d, *r, *z;
>  	int e2=0, e, i, j, l;
> @@ -266,7 +266,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_slots - 1;
>  
>  	do {
>  		*z = y;
> -- 
> 2.21.0
> 

While nothing should change in the case of negative exponents, I also
did an empirical check with DBL_TRUE_MIN to check the space for
denormals, and we have 5 slots to spare.

In the above expression for max_exp_slots, max_exp alone does not
cover the full number of times a digit could be added by halving when
applying the exponent, due to both denormals extending the range, and
the y *= 0x1p28 above (which ensures we pull off as many digits as
possible in the first loop iteration).

However, for denormals, every extra halving that can happen via the
frexp-normalized exponent is one that can't happen in the mantissa. In
other words, all denormals end in the exact same decimal place as the
smallest normal. So they don't contribute anything to the needed
space.

The multiplication by 0x1p28 above is handled by the +28 in the
expression for max_exp_slots.

And the +8 ensures we round up when dividing by 9 decimal digits per
slot.

Re: [musl] Wrong formatting of doubles?

[Markus Wichmann]

Am Wed, Jul 02, 2025 at 10:58:23AM -0400 schrieb Rich Felker:
> On Tue, Jul 01, 2025 at 09:46:10PM -0400, Rich Felker wrote:
> >  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 max_mant_slots = (max_mant_dig+28)/29 + 1;
> > +	int max_exp_slots = (max_exp+max_mant_dig+28+8)/9;
> > +	int bufsize = max_mant_slots + max_exp_slots;
> >  	uint32_t big[bufsize];
> >  	uint32_t *a, *d, *r, *z;
> >  	int e2=0, e, i, j, l;
> > @@ -266,7 +266,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_slots - 1;
> >  
> >  	do {
> >  		*z = y;
> > -- 
> > 2.21.0
> > 
> 
> While nothing should change in the case of negative exponents, I also
> did an empirical check with DBL_TRUE_MIN to check the space for
> denormals, and we have 5 slots to spare.

That tracks with my calculations. I would just look at it from the
result: The buffer only needs to be big enough to contain the largest
end result, right? The largest decimal expansion occurs at TRUE_MIN,
which is 2^(MIN_EXP - MANT_DIG), which has MANT_DIG - MIN_EXP decimal
places. So you need enough slots for that many digits, +1 for the
pre-radix slot you always have, so in the end

bufsize = (mant_dig - min_exp + 8)/9 + 1

And that results in a little bit less than what you calculate.

type | yours | mine
-----|-------|-----
ld64 | 126   | 121
ld80 | 1835  | 1828
ld128| 1842  | 1833

> 
> The multiplication by 0x1p28 above is handled by the +28 in the
> expression for max_exp_slots.
> 

See, that is an intermediate result. The multiplication is not big
enough to make the mantissa expansion require more pre-radix slots, so
it doesn't change how much space is needed. It only changes how many
slots the mantissa expansion might take, which is not important to
printing TRUE_MIN. It is important to printing numbers with positive
exponents, but only for keeping the right distance from the end of the
buffer.

Ciao,
Markus

Re: Wrong formatting of doubles?

[Michal Terepeta]

[-- Attachment #1: Type: text/plain, Size: 528 bytes --]

Hi Rich,

Did you merge this? Is there anything that I can do to help?

Thanks!

- Michal

On Wed, Jul 2, 2025 at 9:56 AM Michal Terepeta <michalt@google.com> wrote:
>
> On Wed, Jul 2, 2025 at 3:46 AM Rich Felker <dalias@libc.org> wrote:
> > [...]
> >
> > Proposed commit attached. Comments welcome. I'll plan to push tomorrow
> > if there are no objections.
> >
>
> I just tested the patch on our end (64-bit RISC-V) and everything looks good!
> Thank you for a quick response and fix! :D
>
> - Michal

[-- Attachment #2: S/MIME Cryptographic Signature --]
[-- Type: application/pkcs7-signature, Size: 5279 bytes --]

Re: [musl] Wrong formatting of doubles?

[Rich Felker]

On Tue, Jul 08, 2025 at 09:43:45AM +0200, Michal Terepeta wrote:
> Hi Rich,
> 
> Did you merge this? Is there anything that I can do to help?

No worries, I just realized I have the arm __getauxval patch discussed
in the SME thread also queued and was waiting to see what to do with
that and whether I should rebase without it before pushing. But I
haven't forgotten it, and it should be fine to use the version I
posted. If the arm thing isn't resolved quickly I'll just rebase and
push without it for now.

Rich


> Thanks!
> 
> - Michal
> 
> On Wed, Jul 2, 2025 at 9:56 AM Michal Terepeta <michalt@google.com> wrote:
> >
> > On Wed, Jul 2, 2025 at 3:46 AM Rich Felker <dalias@libc.org> wrote:
> > > [...]
> > >
> > > Proposed commit attached. Comments welcome. I'll plan to push tomorrow
> > > if there are no objections.
> > >
> >
> > I just tested the patch on our end (64-bit RISC-V) and everything looks good!
> > Thank you for a quick response and fix! :D
> >
> > - Michal