Re: Wrong formatting of doubles?

[Rich Felker <dalias@libc.org>]

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