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

[Peter Ammon <corydoras@ridiculousfish.com>]

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