[musl] [PATCH] Fix printf hex float formatting with precision

[musl] [PATCH] Fix printf hex float formatting with precision

[Peter Ammon]

printf hex formatting with precision may emit excess digits on targets where long double is an alias of double. For example, on ARMv7, the expression `printf("%.12a", M_PI)` outputs 13 digits past the decimal, instead of 12.

I believe the cause is a bogus rounding calculation in fmt_fp:

if (p<0 || p>=LDBL_MANT_DIG/4-1) re=0;
else re=LDBL_MANT_DIG/4-1-p;
if (re) {
    round *= 1<<(LDBL_MANT_DIG%4);
    while (re--) round*=16;
    ...

I wasn't able to justify the calculation of `re`; I think it suffers from trying to round in terms of number of hex digits, which is tricky because the number of fractional mantissa bits may not be a multiple of 4.

I propose to fix it by working in bits instead of hex digits, as follows:

if (p>=0 && p<(LDBL_MANT_DIG-1+3)/4) {
    int re = LDBL_MANT_DIG-1-(p*4);
    long double round = 1ULL<<re;
...

This is justified as follows:

- The value to format has been scaled to the range [1, 2), so there is exactly one bit before the radix. Thus, the fractional portion of the mantissa may be printed at full precision with LDBL_MANT_DIG-1, divided by 4, rounding up; thus `(LDBL_MANT_DIG-1+3)/4`
- The caller has requested a precision of p hex digits, so p*4 bits.
- `re` is then the number of bits to round off, as LDBL_MANT_DIG-1-(p*4)
- Thus we compute `round` as 2**re. Adding `round` will lose `re` bits of precision, rounding per the fp mode; subtracting this again will round the original value.

I also removed an initial factor of 8 from the `round` as it seemed incorrect to me; for example we may want to round only the last bit, in which case the min round of 8 would be too big.

I am by no means a numerics expert, so I very much welcome another pair of eyes. I will be happy to contribute a test.

---
 src/stdio/vfprintf.c | 12 +++---------
 1 file changed, 3 insertions(+), 9 deletions(-)

diff --git a/src/stdio/vfprintf.c b/src/stdio/vfprintf.c
index 497c5e19..88ecbfa2 100644
--- a/src/stdio/vfprintf.c
+++ b/src/stdio/vfprintf.c
@@ -211,18 +211,12 @@ static int fmt_fp(FILE *f, long double y, int w, int p, int fl, int t)
 	if (y) e2--;

 	if ((t|32)=='a') {
-		long double round = 8.0;
-		int re;
-
 		if (t&32) prefix += 9;
 		pl += 2;

-		if (p<0 || p>=LDBL_MANT_DIG/4-1) re=0;
-		else re=LDBL_MANT_DIG/4-1-p;
-
-		if (re) {
-			round *= 1<<(LDBL_MANT_DIG%4);
-			while (re--) round*=16;
+		if (p>=0 && p<(LDBL_MANT_DIG-1+3)/4) {
+			int re = LDBL_MANT_DIG-1-(p*4);
+			long double round = 1ULL<<re;
 			if (*prefix=='-') {
 				y=-y;
 				y-=round;
--
2.39.2

Re: [musl] [PATCH] Fix printf hex float formatting with precision

[Rich Felker]

On Thu, Apr 11, 2024 at 06:17:25PM -0700, Peter Ammon wrote:
> printf hex formatting with precision may emit excess digits on
> targets where long double is an alias of double. For example, on
> ARMv7, the expression `printf("%.12a", M_PI)` outputs 13 digits past
> the decimal, instead of 12.
> 
> I believe the cause is a bogus rounding calculation in fmt_fp:
> 
> if (p<0 || p>=LDBL_MANT_DIG/4-1) re=0;
> else re=LDBL_MANT_DIG/4-1-p;
> if (re) {
>     round *= 1<<(LDBL_MANT_DIG%4);
>     while (re--) round*=16;
>     ...
> 
> I wasn't able to justify the calculation of `re`; I think it suffers
> from trying to round in terms of number of hex digits, which is
> tricky because the number of fractional mantissa bits may not be a
> multiple of 4.
> 
> I propose to fix it by working in bits instead of hex digits, as follows:
> 
> if (p>=0 && p<(LDBL_MANT_DIG-1+3)/4) {
>     int re = LDBL_MANT_DIG-1-(p*4);
>     long double round = 1ULL<<re;
> ....
> 
> This is justified as follows:
> 
> - The value to format has been scaled to the range [1, 2), so there
> is exactly one bit before the radix. Thus, the fractional portion of
> the mantissa may be printed at full precision with LDBL_MANT_DIG-1,
> divided by 4, rounding up; thus `(LDBL_MANT_DIG-1+3)/4`
> - The caller has requested a precision of p hex digits, so p*4 bits.
> - `re` is then the number of bits to round off, as LDBL_MANT_DIG-1-(p*4)
> - Thus we compute `round` as 2**re. Adding `round` will lose `re`
> bits of precision, rounding per the fp mode; subtracting this again
> will round the original value.
> 
> I also removed an initial factor of 8 from the `round` as it seemed
> incorrect to me; for example we may want to round only the last bit,
> in which case the min round of 8 would be too big.
> 
> I am by no means a numerics expert, so I very much welcome another
> pair of eyes. I will be happy to contribute a test.

Aside from one minor fix, this code dates back to the initial check-in
of musl, and based on the behavior and the initial value 8 for round,
I'm pretty sure it was written assuming we want the digit before the
radix point to be in the range [8,F] rather than 1, to pack the max
amount of precision into the output. For some reason, I decided this
was not the best behavior, but must have failed to compensate for that
in the rounding logic. There may be pre-0.5.0 trees somewhere I could
excavate to figure out exactly what happened, but that's my best
guess.

FWIW the bug seems to be fixed by making it do y*=8; e2-=3;

Offhand without delving deep into it, your fix sounds correct. It
could probably be factored as two parts, first just fixing the bug
(I think the only actual bug is the 8), and a second switching to your
better algorithm for computing the rounding term.

Rich

Re: [musl] [PATCH] Fix printf hex float formatting with precision

[Rich Felker]

On Thu, Apr 11, 2024 at 06:17:25PM -0700, Peter Ammon wrote:
> if (p>=0 && p<(LDBL_MANT_DIG-1+3)/4) {
>     int re = LDBL_MANT_DIG-1-(p*4);
>     long double round = 1ULL<<re;

This expression overflows. re can be as large as 108 but 1ULL has
64-bit type.

There's probably some reasonable way to write it that works for
reasonable sizes of long double (I think it's safe to assume IEEE quad
is the max that will ever be supported), but I presume I wrote it with
the original inefficient loop to be fully general to arbitrary
precision.

It could just be written to use scalbn, I think. At the time I
probably was trying to avoid dependency on a libm that could have been
separate from libc (bad historical thing). IIRC there was a time when
frexp was not used either.

Rich

Re: [musl] [PATCH] Fix printf hex float formatting with precision

[Rich Felker]

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

On Fri, Apr 12, 2024 at 03:57:28PM -0400, Rich Felker wrote:
> On Thu, Apr 11, 2024 at 06:17:25PM -0700, Peter Ammon wrote:
> > if (p>=0 && p<(LDBL_MANT_DIG-1+3)/4) {
> >     int re = LDBL_MANT_DIG-1-(p*4);
> >     long double round = 1ULL<<re;
> 
> This expression overflows. re can be as large as 108 but 1ULL has
> 64-bit type.
> 
> There's probably some reasonable way to write it that works for
> reasonable sizes of long double (I think it's safe to assume IEEE quad
> is the max that will ever be supported), but I presume I wrote it with
> the original inefficient loop to be fully general to arbitrary
> precision.
> 
> It could just be written to use scalbn, I think. At the time I
> probably was trying to avoid dependency on a libm that could have been
> separate from libc (bad historical thing). IIRC there was a time when
> frexp was not used either.

Does the attached look ok?

Rich

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

diff --git a/src/stdio/vfprintf.c b/src/stdio/vfprintf.c
index 497c5e19..dc648e7e 100644
--- a/src/stdio/vfprintf.c
+++ b/src/stdio/vfprintf.c
@@ -211,18 +211,11 @@ static int fmt_fp(FILE *f, long double y, int w, int p, int fl, int t)
 	if (y) e2--;
 
 	if ((t|32)=='a') {
-		long double round = 8.0;
-		int re;
-
 		if (t&32) prefix += 9;
 		pl += 2;
 
-		if (p<0 || p>=LDBL_MANT_DIG/4-1) re=0;
-		else re=LDBL_MANT_DIG/4-1-p;
-
-		if (re) {
-			round *= 1<<(LDBL_MANT_DIG%4);
-			while (re--) round*=16;
+		if (p>=0 && p<(LDBL_MANT_DIG-1+3)/4) {
+			double round = scalbn(1, LDBL_MANT_DIG-1-(p*4));
 			if (*prefix=='-') {
 				y=-y;
 				y-=round;

Re: [musl] [PATCH] Fix printf hex float formatting with precision

[Peter Ammon]

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

Good catch on the overflow. Yes, that patch looks correct to me and I
confirmed it fixes the issue on ARMv7.

> On Apr 12, 2024, at 4:33 PM, Rich Felker <dalias@libc.org> wrote:
> 
> On Fri, Apr 12, 2024 at 03:57:28PM -0400, Rich Felker wrote:
>> On Thu, Apr 11, 2024 at 06:17:25PM -0700, Peter Ammon wrote:
>>> if (p>=0 && p<(LDBL_MANT_DIG-1+3)/4) {
>>>    int re = LDBL_MANT_DIG-1-(p*4);
>>>    long double round = 1ULL<<re;
>> 
>> This expression overflows. re can be as large as 108 but 1ULL has
>> 64-bit type.
>> 
>> There's probably some reasonable way to write it that works for
>> reasonable sizes of long double (I think it's safe to assume IEEE quad
>> is the max that will ever be supported), but I presume I wrote it with
>> the original inefficient loop to be fully general to arbitrary
>> precision.
>> 
>> It could just be written to use scalbn, I think. At the time I
>> probably was trying to avoid dependency on a libm that could have been
>> separate from libc (bad historical thing). IIRC there was a time when
>> frexp was not used either.
> 
> Does the attached look ok?
> 
> Rich
> <pct_a.diff>


[-- Attachment #2: Type: text/html, Size: 4727 bytes --]