printf issues

printf issues

[Morten Welinder]

I am seeing differences in how printf("%.200Lg",val) works between musl code
and glibc.

Here are some samples.  The top line is how musl prints, the bottom is glibc.

-65878995336522048.000000000
-65878995336522048

1954675876964773.500000000
1954675876964773.5

3953605802361882.000000000
3953605802361882


Unrelatedly, from function fmt_fp:

#define CONCAT2(x,y) x ## y
#define CONCAT(x,y) CONCAT2(x,y)
[...]
            long double round = CONCAT(0x1p,LDBL_MANT_DIG);

That code is cute as a Hello Kitty door knocker, but really?  Let's hope nobody
gets the urge to define LDBL_MANT_DIG as 0100 or (80-16) or some such.
The first case will still compile, but get the wrong result.

Morten

Re: printf issues

[Rich Felker]

On Fri, Apr 04, 2014 at 09:32:18AM -0400, Morten Welinder wrote:
> I am seeing differences in how printf("%.200Lg",val) works between musl code
> and glibc.
> 
> Here are some samples.  The top line is how musl prints, the bottom is glibc.
> 
> -65878995336522048.000000000
> -65878995336522048
> 
> 1954675876964773.500000000
> 1954675876964773.5
> 
> 3953605802361882.000000000
> 3953605802361882

Confirmed that this bug exists. It seems to be an error in how the
trailing zero removal logic works in certain cases, but I didn't
figure out the pattern immediately except that it seems to be
associated with large-magnitude numbers like yours. I'll run some
better tests in a little bit.

> Unrelatedly, from function fmt_fp:
> 
> #define CONCAT2(x,y) x ## y
> #define CONCAT(x,y) CONCAT2(x,y)
> [...]
>             long double round = CONCAT(0x1p,LDBL_MANT_DIG);
> 
> That code is cute as a Hello Kitty door knocker, but really?  Let's hope nobody
> gets the urge to define LDBL_MANT_DIG as 0100 or (80-16) or some such.
> The first case will still compile, but get the wrong result.

Yes, this assumes our headers which define it as an unadorned decimal
constant. Do you have any ideas for a clean way to avoid this
assumption without having to compute the value at runtime?

Rich

Re: printf issues

[Szabolcs Nagy]

* Morten Welinder <mwelinder@gmail.com> [2014-04-04 09:32:18 -0400]:
> Unrelatedly, from function fmt_fp:
> 
> #define CONCAT2(x,y) x ## y
> #define CONCAT(x,y) CONCAT2(x,y)
> [...]
>             long double round = CONCAT(0x1p,LDBL_MANT_DIG);
> 
> That code is cute as a Hello Kitty door knocker, but really?  Let's hope nobody
> gets the urge to define LDBL_MANT_DIG as 0100 or (80-16) or some such.

before you can mock libc code you need to educate yourself

LDBL_MANT_DIG is an identifier (macro definition) that is
defined by the libc itself (in arch/$ARCH/bits/float.h),
not by the application, nor by the compiler

the libc can internally rely on a different contract about
the identifiers it defines, than application code that has
to rely on the external contracts specified by ISO C

so outside the implementation you shouldn't use such a code
because there is no guarantee how LDBL_MANT_DIG is defined,
but internally the libc knows how it defined its symbols and
thus the code is perfectly reasonable

(this distinction between public interface contracts and
implementation internal interface contracts is often
misunderstood: so don't copy-paste code between libc and
application code without thinking and don't try to reason
about APIs/ABIs used inside the libc based on the public
specification of those APIs/ABIs)

Re: printf issues

[Morten Welinder]

> before you can mock libc code you need to educate yourself

I did and that's why I called the code "cute", not "wrong".  But if you read
the porting documentation

    http://brightrain.aerifal.cx/~niklata/PORTING
    http://www.openwall.com/lists/musl/2012/07/08/1

you will notice that nowhere does it warn that defining LDBL_MANT_DIG
as anything but a base-10 constant may cause printf-rounding to fail.

> Do you have any ideas for a clean way to avoid this
> assumption without having to compute the value at runtime?

I don't know if ldexpl will get constant folded by the compiler, but if not,
I think (2.0L/LDBL_EPSILON) ought to work as a replacement.  It's not
as likely to get prices at the obfuscated C contents, though.

Morten



On Fri, Apr 4, 2014 at 10:15 AM, Szabolcs Nagy <nsz@port70.net> wrote:
> * Morten Welinder <mwelinder@gmail.com> [2014-04-04 09:32:18 -0400]:
>> Unrelatedly, from function fmt_fp:
>>
>> #define CONCAT2(x,y) x ## y
>> #define CONCAT(x,y) CONCAT2(x,y)
>> [...]
>>             long double round = CONCAT(0x1p,LDBL_MANT_DIG);
>>
>> That code is cute as a Hello Kitty door knocker, but really?  Let's hope nobody
>> gets the urge to define LDBL_MANT_DIG as 0100 or (80-16) or some such.
>
> before you can mock libc code you need to educate yourself
>
> LDBL_MANT_DIG is an identifier (macro definition) that is
> defined by the libc itself (in arch/$ARCH/bits/float.h),
> not by the application, nor by the compiler
>
> the libc can internally rely on a different contract about
> the identifiers it defines, than application code that has
> to rely on the external contracts specified by ISO C
>
> so outside the implementation you shouldn't use such a code
> because there is no guarantee how LDBL_MANT_DIG is defined,
> but internally the libc knows how it defined its symbols and
> thus the code is perfectly reasonable
>
> (this distinction between public interface contracts and
> implementation internal interface contracts is often
> misunderstood: so don't copy-paste code between libc and
> application code without thinking and don't try to reason
> about APIs/ABIs used inside the libc based on the public
> specification of those APIs/ABIs)

Re: printf issues

[Szabolcs Nagy]

* Morten Welinder <mwelinder@gmail.com> [2014-04-04 10:35:00 -0400]:
> > Do you have any ideas for a clean way to avoid this
> > assumption without having to compute the value at runtime?
> 
> I don't know if ldexpl will get constant folded by the compiler, but if not,
> I think (2.0L/LDBL_EPSILON) ought to work as a replacement.  It's not
> as likely to get prices at the obfuscated C contents, though.

ldexpl is not ok
(cannot get constant folded in free-standing mode unless
we use __builtin_ldexpl on compilers that support that)

the epsilon based solution looks good though

Re: printf issues

[Rich Felker]

On Fri, Apr 04, 2014 at 10:35:00AM -0400, Morten Welinder wrote:
> > before you can mock libc code you need to educate yourself
> 
> I did and that's why I called the code "cute", not "wrong".  But if you read
> the porting documentation
> 
>     http://brightrain.aerifal.cx/~niklata/PORTING
>     http://www.openwall.com/lists/musl/2012/07/08/1
> 
> you will notice that nowhere does it warn that defining LDBL_MANT_DIG
> as anything but a base-10 constant may cause printf-rounding to fail.

Good point.

> > Do you have any ideas for a clean way to avoid this
> > assumption without having to compute the value at runtime?
> 
> I don't know if ldexpl will get constant folded by the compiler, but if not,
> I think (2.0L/LDBL_EPSILON) ought to work as a replacement.  It's not
> as likely to get prices at the obfuscated C contents, though.

Thanks, I think that's exactly the right solution. FWIW, I _would_
like this code to be easily adaptable for use outside libc if somebody
wants it, so eliminating implementation-internal assumptions like this
is nice.

Rich

Re: printf issues

[Morten Welinder]

> I _would_ like this code to be easily adaptable for use outside libc
> if somebody wants it

FYI, I have been doing just that for Gnumeric in a variant that always
rounds ties away from zero.  Two changes would help with making
the code fit seamlessly into other environments.

1. Make "i" in fmt_fp unsigned.  It's used in connection with
    unsigned values only.

2. Make "char *s" used to hold "NAN" etc. "const char *s".

Neither of these should make any difference in what the function
actually does.

I have run tens of millions random numbers through this function
looking for differences between it and glibc.  The extra 0s from "%g"
is the only problem observed.


It looks like the LDBL_EPSILON version could be used in

    roundl.c
    modfl.c
    ceill.c
    floorl.c

in the definition of TOINT instead of enumerating choices for
LDBL_MANT_DIG.  It's basically the same thing going on
there.

While I was looking for that, I noticed that this modfl fallback looks
problematic.  Even if long double and double are the same thing
under the hood, I don't think you can cast pointers like that and
assume it works.  It needs a temporary.

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double modfl(long double x, long double *iptr)
{
return modf(x, (double *)iptr);
}



Morten





On Fri, Apr 4, 2014 at 11:07 AM, Rich Felker <dalias@aerifal.cx> wrote:
> On Fri, Apr 04, 2014 at 10:35:00AM -0400, Morten Welinder wrote:
>> > before you can mock libc code you need to educate yourself
>>
>> I did and that's why I called the code "cute", not "wrong".  But if you read
>> the porting documentation
>>
>>     http://brightrain.aerifal.cx/~niklata/PORTING
>>     http://www.openwall.com/lists/musl/2012/07/08/1
>>
>> you will notice that nowhere does it warn that defining LDBL_MANT_DIG
>> as anything but a base-10 constant may cause printf-rounding to fail.
>
> Good point.
>
>> > Do you have any ideas for a clean way to avoid this
>> > assumption without having to compute the value at runtime?
>>
>> I don't know if ldexpl will get constant folded by the compiler, but if not,
>> I think (2.0L/LDBL_EPSILON) ought to work as a replacement.  It's not
>> as likely to get prices at the obfuscated C contents, though.
>
> Thanks, I think that's exactly the right solution. FWIW, I _would_
> like this code to be easily adaptable for use outside libc if somebody
> wants it, so eliminating implementation-internal assumptions like this
> is nice.
>
> Rich

Re: printf issues

[Szabolcs Nagy]

* Morten Welinder <mwelinder@gmail.com> [2014-04-04 13:42:30 -0400]:
> It looks like the LDBL_EPSILON version could be used in
> 
>     roundl.c
>     modfl.c
>     ceill.c
>     floorl.c
> 
> in the definition of TOINT instead of enumerating choices for
> LDBL_MANT_DIG.  It's basically the same thing going on

yes, that would be a bit nicer
(although other long double formats won't be supported anytime soon)

(note that in the future these implementations may need to change
the current versions raise inexact flag if result!=input, but the
next version of the floating-point extension standard for c
will require suppressing inexact, which i dont know how to do
with simple arithmetics efficiently without accessing the fenv..)

> While I was looking for that, I noticed that this modfl fallback looks
> problematic.  Even if long double and double are the same thing
> under the hood, I don't think you can cast pointers like that and
> assume it works.  It needs a temporary.
> 
> #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
> long double modfl(long double x, long double *iptr)
> {
> return modf(x, (double *)iptr);
> }

yes, this is an aliasing violation, nice catch

the original idea was to allow tail call opt for these wrappers,
so they are a single branch instruction, we should fix it but
i think we can rely on that the ptr representations are the same:

long double modfl(long double x, long double *iptr)
{
	union {long double *ld; double *d;} u = {iptr};
	return modf(x, u.d);
}

Re: printf issues

[Morten Welinder]

In fmt_fmt, the rounding decision is done using this test:

            /* Decide whether to round by probing round+small */
            if (round+small != round) { ...

Why is this done with long double?

The reason I ask is that the Valgrind situation improves a lot if
this is done with doubles.

(Valgrind situation: Valgrind emulates long doubles, poorly, by using
simple doubles.  See, for example, https://bugs.kde.org/show_bug.cgi?id=164298)

Morten




On Fri, Apr 4, 2014 at 2:54 PM, Szabolcs Nagy <nsz@port70.net> wrote:
> * Morten Welinder <mwelinder@gmail.com> [2014-04-04 13:42:30 -0400]:
>> It looks like the LDBL_EPSILON version could be used in
>>
>>     roundl.c
>>     modfl.c
>>     ceill.c
>>     floorl.c
>>
>> in the definition of TOINT instead of enumerating choices for
>> LDBL_MANT_DIG.  It's basically the same thing going on
>
> yes, that would be a bit nicer
> (although other long double formats won't be supported anytime soon)
>
> (note that in the future these implementations may need to change
> the current versions raise inexact flag if result!=input, but the
> next version of the floating-point extension standard for c
> will require suppressing inexact, which i dont know how to do
> with simple arithmetics efficiently without accessing the fenv..)
>
>> While I was looking for that, I noticed that this modfl fallback looks
>> problematic.  Even if long double and double are the same thing
>> under the hood, I don't think you can cast pointers like that and
>> assume it works.  It needs a temporary.
>>
>> #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
>> long double modfl(long double x, long double *iptr)
>> {
>> return modf(x, (double *)iptr);
>> }
>
> yes, this is an aliasing violation, nice catch
>
> the original idea was to allow tail call opt for these wrappers,
> so they are a single branch instruction, we should fix it but
> i think we can rely on that the ptr representations are the same:
>
> long double modfl(long double x, long double *iptr)
> {
>         union {long double *ld; double *d;} u = {iptr};
>         return modf(x, u.d);
> }

Re: printf issues

[Morten Welinder]

Another printf issue has shown up, this time with memory corruption.

    printf ("%.3E\n", 999999999.0);

The rounding test correctly decides that it needs to round this value
up to 1E+09.  It is, however, utterly unprepared for having nowhere to
put the carry.  It happily accesses and changes one or more elements
before the one that held 999999999.

Morten




On Fri, Apr 4, 2014 at 4:01 PM, Morten Welinder <mwelinder@gmail.com> wrote:
> In fmt_fmt, the rounding decision is done using this test:
>
>             /* Decide whether to round by probing round+small */
>             if (round+small != round) { ...
>
> Why is this done with long double?
>
> The reason I ask is that the Valgrind situation improves a lot if
> this is done with doubles.
>
> (Valgrind situation: Valgrind emulates long doubles, poorly, by using
> simple doubles.  See, for example, https://bugs.kde.org/show_bug.cgi?id=164298)
>
> Morten
>
>
>
>
> On Fri, Apr 4, 2014 at 2:54 PM, Szabolcs Nagy <nsz@port70.net> wrote:
>> * Morten Welinder <mwelinder@gmail.com> [2014-04-04 13:42:30 -0400]:
>>> It looks like the LDBL_EPSILON version could be used in
>>>
>>>     roundl.c
>>>     modfl.c
>>>     ceill.c
>>>     floorl.c
>>>
>>> in the definition of TOINT instead of enumerating choices for
>>> LDBL_MANT_DIG.  It's basically the same thing going on
>>
>> yes, that would be a bit nicer
>> (although other long double formats won't be supported anytime soon)
>>
>> (note that in the future these implementations may need to change
>> the current versions raise inexact flag if result!=input, but the
>> next version of the floating-point extension standard for c
>> will require suppressing inexact, which i dont know how to do
>> with simple arithmetics efficiently without accessing the fenv..)
>>
>>> While I was looking for that, I noticed that this modfl fallback looks
>>> problematic.  Even if long double and double are the same thing
>>> under the hood, I don't think you can cast pointers like that and
>>> assume it works.  It needs a temporary.
>>>
>>> #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
>>> long double modfl(long double x, long double *iptr)
>>> {
>>> return modf(x, (double *)iptr);
>>> }
>>
>> yes, this is an aliasing violation, nice catch
>>
>> the original idea was to allow tail call opt for these wrappers,
>> so they are a single branch instruction, we should fix it but
>> i think we can rely on that the ptr representations are the same:
>>
>> long double modfl(long double x, long double *iptr)
>> {
>>         union {long double *ld; double *d;} u = {iptr};
>>         return modf(x, u.d);
>> }

Re: printf issues

[Rich Felker]

On Fri, Apr 04, 2014 at 04:22:46PM -0400, Morten Welinder wrote:
> Another printf issue has shown up, this time with memory corruption.
> 
>     printf ("%.3E\n", 999999999.0);
> 
> The rounding test correctly decides that it needs to round this value
> up to 1E+09.  It is, however, utterly unprepared for having nowhere to
> put the carry.  It happily accesses and changes one or more elements
> before the one that held 999999999.

I suspect this may be true; if so, it's a very nice catch. Were you
able to determine what data it clobbers (in practice; obviously this
is compiler-specific) and whether the clobber has any observable
effects?

Rich

Re: printf issues

[Morten Welinder]

> Were you able to determine what data it clobbers (in practice;
> obviously this is compiler-specific) and whether the clobber
> has any observable effects?

It clobbers uninitialized parts of "big".  If you add

    for (i = 0; i < sizeof(big)/sizeof(big[0]); i++) big[i] = 12345678;

then it will consistently print "1.23E+16" which is a bit off, :-)  If
you instead
initialize like this:

    for (i = 0; i < sizeof(big)/sizeof(big[0]); i++) big[i] = 999999999;

then I get "1E+15939" which is fairly impressive.  Also, in this case it will
clobber whatever happened to come before "big".

Morten



On Fri, Apr 4, 2014 at 5:08 PM, Rich Felker <dalias@aerifal.cx> wrote:
> On Fri, Apr 04, 2014 at 04:22:46PM -0400, Morten Welinder wrote:
>> Another printf issue has shown up, this time with memory corruption.
>>
>>     printf ("%.3E\n", 999999999.0);
>>
>> The rounding test correctly decides that it needs to round this value
>> up to 1E+09.  It is, however, utterly unprepared for having nowhere to
>> put the carry.  It happily accesses and changes one or more elements
>> before the one that held 999999999.
>
> I suspect this may be true; if so, it's a very nice catch. Were you
> able to determine what data it clobbers (in practice; obviously this
> is compiler-specific) and whether the clobber has any observable
> effects?
>
> Rich

Re: printf issues

[Morten Welinder]

I *think* the right fix is to add the following "if' statement into
the rounding loop:

                while (*d > 999999999) {
                    *d--=0;
                    if (d < a) *--a = 0;
                    (*d)++;
                }

This also ought to make the d<a test afterwards unnecessary.  But
more tests would be better.

M.




On Fri, Apr 4, 2014 at 6:50 PM, Morten Welinder <mwelinder@gmail.com> wrote:
>> Were you able to determine what data it clobbers (in practice;
>> obviously this is compiler-specific) and whether the clobber
>> has any observable effects?
>
> It clobbers uninitialized parts of "big".  If you add
>
>     for (i = 0; i < sizeof(big)/sizeof(big[0]); i++) big[i] = 12345678;
>
> then it will consistently print "1.23E+16" which is a bit off, :-)  If
> you instead
> initialize like this:
>
>     for (i = 0; i < sizeof(big)/sizeof(big[0]); i++) big[i] = 999999999;
>
> then I get "1E+15939" which is fairly impressive.  Also, in this case it will
> clobber whatever happened to come before "big".
>
> Morten
>
>
>
> On Fri, Apr 4, 2014 at 5:08 PM, Rich Felker <dalias@aerifal.cx> wrote:
>> On Fri, Apr 04, 2014 at 04:22:46PM -0400, Morten Welinder wrote:
>>> Another printf issue has shown up, this time with memory corruption.
>>>
>>>     printf ("%.3E\n", 999999999.0);
>>>
>>> The rounding test correctly decides that it needs to round this value
>>> up to 1E+09.  It is, however, utterly unprepared for having nowhere to
>>> put the carry.  It happily accesses and changes one or more elements
>>> before the one that held 999999999.
>>
>> I suspect this may be true; if so, it's a very nice catch. Were you
>> able to determine what data it clobbers (in practice; obviously this
>> is compiler-specific) and whether the clobber has any observable
>> effects?
>>
>> Rich

Re: printf issues

[Rich Felker]

On Fri, Apr 04, 2014 at 08:01:00PM -0400, Morten Welinder wrote:
> I *think* the right fix is to add the following "if' statement into
> the rounding loop:
> 
>                 while (*d > 999999999) {
>                     *d--=0;
>                     if (d < a) *--a = 0;
>                     (*d)++;
>                 }
> 
> This also ought to make the d<a test afterwards unnecessary.  But
> more tests would be better.

This looks right.

Are you sure a==big is not possible at this point? I *think* it's
impossible, since a==big can only happen if e2<0, in which case
applying the exponent is going to reduce the first slot so that
rounding could never cause it to overflow. Do you agree with this
analysis?

Rich

Re: printf issues

[Rich Felker]

On Fri, Apr 04, 2014 at 08:01:00PM -0400, Morten Welinder wrote:
> I *think* the right fix is to add the following "if' statement into
> the rounding loop:
> 
>                 while (*d > 999999999) {
>                     *d--=0;
>                     if (d < a) *--a = 0;
>                     (*d)++;
>                 }
> 
> This also ought to make the d<a test afterwards unnecessary.  But
> more tests would be better.

After re-studying the code, I think this is the correct fix, and I've
committed the fix to git. I also fixed the bug where %g failed to trim
trailing zeros.

Rich

Re: printf issues

[Morten Welinder]

I am seeing another case of printing the wrong result.  (This is with
the pre-fix code, but I don't think that matters.)

For 0xc.301316272b908p+49L (aka 6861116509411105.0L) with format
"%.15Lg" I get

musl: 6.86111650941111e+15
glibc: 6.8611165094111e+15

Musl has rounded to odd here for a midpoint case.

Morten







On Mon, Apr 7, 2014 at 3:29 AM, Rich Felker <dalias@aerifal.cx> wrote:
> On Fri, Apr 04, 2014 at 08:01:00PM -0400, Morten Welinder wrote:
>> I *think* the right fix is to add the following "if' statement into
>> the rounding loop:
>>
>>                 while (*d > 999999999) {
>>                     *d--=0;
>>                     if (d < a) *--a = 0;
>>                     (*d)++;
>>                 }
>>
>> This also ought to make the d<a test afterwards unnecessary.  But
>> more tests would be better.
>
> After re-studying the code, I think this is the correct fix, and I've
> committed the fix to git. I also fixed the bug where %g failed to trim
> trailing zeros.
>
> Rich

Re: printf issues

[Morten Welinder]

The frequency of this problem is something like 1 in 5e6.
Observations:

*  I only seem to be able to trigger it for %g even though all my
   samples print in "e" form.

* The numbers are all roughly the same size: 1e15

* The numbers are all integers ending in 05.  (Except two cases
   where the 05 is followed by zeros.)

* The precision is always just below the value that would have
  make an exact representation.

Morten



Test #1165601 at precision 15:
-4.15655192121011e+15
-4.1565519212101e+15
-0xe.c45ca8c112f9p+48
-4156551921210105

Test #7022240 at precision 15:
-8.18656081955811e+15
-8.1865608195581e+15
-0xe.8ad11cca7b6c8p+49
-8186560819558105

Test #7718112 at precision 15:
-3.57044516323631e+15
-3.5704451632363e+15
-0xc.af4d0ba50bd1p+48
-3570445163236305

Test #7772566 at precision 15:
2.50163300941741e+15
2.5016330094174e+15
0x8.e338d2e27cbdp+48
2501633009417405

Test #8135560 at precision 15:
-8.75451694516931e+15
-8.7545169451693e+15
-0xf.8d17e85935cc8p+49
-8754516945169305

Test #11526318 at precision 15:
7.64406482711661e+15
7.6440648271166e+15
0xd.941e4454861e8p+49
7644064827116605

Test #13957100 at precision 14:
-9.7205430380441e+14
-9.720543038044e+14
-0xd.d0501dec2fd4p+46
-972054303804405

Test #20008327 at precision 15:
3.37380293784171e+15
3.3738029378417e+15
0xb.fc74b5a6f829p+48
3373802937841705

Test #23073645 at precision 14:
-5.1833238325691e+15
-5.183323832569e+15
-0x9.351a4fe5e66dp+49
-5183323832569050

Test #25660410 at precision 15:
3.96328205296061e+15
3.9632820529606e+15
0xe.149582e9515dp+48
3963282052960605

Test #29106716 at precision 13:
8.042000473751e+15
8.04200047375e+15
0xe.491412c2bdf2p+49
8042000473750500

Test #35240073 at precision 15:
1.96456678790721e+15
1.9645667879072e+15
0xd.f586b30fbd0ap+47
1964566787907205

Test #35298996 at precision 15:
6.67034256219981e+15
6.6703425621998e+15
0xb.d95213799c7e8p+49
6670342562199805

Test #36728927 at precision 13:
-2.920984132831e+13
-2.92098413283e+13
-0xd.48791bb0588p+41
-29209841328305

Test #40382272 at precision 15:
8.82887106290361e+15
8.8288710629036e+15
0xf.aee7ddbc6d9a8p+49
8828871062903605

Test #44394631 at precision 14:
-2.4546341550581e+14
-2.454634155058e+14
-0xd.f3f67afc38dp+44
-245463415505805

Test #56164412 at precision 15:
-5.26534878693031e+15
-5.2653487869303e+15
-0x9.5a67460821408p+49
-5265348786930305

Test #58309890 at precision 15:
-3.54820922268461e+15
-3.5482092226846e+15
-0xc.9b13d64e8fbdp+48
-3548209222684605

Test #64933676 at precision 15:
5.61427609297841e+15
5.6142760929784e+15
0x9.f913c218b2728p+49
5614276092978405

Test #70560963 at precision 15:
-2.75797952567011e+15
-2.7579795256701e+15
-0x9.cc5e25ece4d9p+48
-2757979525670105

Test #80659572 at precision 15:
1.99612924552141e+15
1.9961292455214e+15
0xe.2ef01d35cbfap+47
1996129245521405

On Mon, Apr 7, 2014 at 9:40 AM, Morten Welinder <mwelinder@gmail.com> wrote:
> I am seeing another case of printing the wrong result.  (This is with
> the pre-fix code, but I don't think that matters.)
>
> For 0xc.301316272b908p+49L (aka 6861116509411105.0L) with format
> "%.15Lg" I get
>
> musl: 6.86111650941111e+15
> glibc: 6.8611165094111e+15
>
> Musl has rounded to odd here for a midpoint case.
>
> Morten
>
>
>
>
>
>
>
> On Mon, Apr 7, 2014 at 3:29 AM, Rich Felker <dalias@aerifal.cx> wrote:
>> On Fri, Apr 04, 2014 at 08:01:00PM -0400, Morten Welinder wrote:
>>> I *think* the right fix is to add the following "if' statement into
>>> the rounding loop:
>>>
>>>                 while (*d > 999999999) {
>>>                     *d--=0;
>>>                     if (d < a) *--a = 0;
>>>                     (*d)++;
>>>                 }
>>>
>>> This also ought to make the d<a test afterwards unnecessary.  But
>>> more tests would be better.
>>
>> After re-studying the code, I think this is the correct fix, and I've
>> committed the fix to git. I also fixed the bug where %g failed to trim
>> trailing zeros.
>>
>> Rich

Re: printf issues

[Szabolcs Nagy]

* Morten Welinder <mwelinder@gmail.com> [2014-04-07 10:13:26 -0400]:
> The frequency of this problem is something like 1 in 5e6.
> Observations:
> 
> *  I only seem to be able to trigger it for %g even though all my
>    samples print in "e" form.
> 
> * The numbers are all roughly the same size: 1e15
> 
> * The numbers are all integers ending in 05.  (Except two cases
>    where the 05 is followed by zeros.)
> 
> * The precision is always just below the value that would have
>   make an exact representation.

printf("%.12g\n", 1000000000005.0);
printf("%.11g\n", 500000000045.0);
printf("%.11g\n", 275000000025.0);

prints

1.00000000001e+12
5.0000000005e+11
2.7500000003e+11

in fmt_fp
	if (x || d+1!=z) {
		long double round = CONCAT(0x1p,LDBL_MANT_DIG);
		long double small;
		if (*d/i & 1) round += 2;
		if (x<i/2) small=0x0.8p0;
		else if (x==i/2 && d+1==z) small=0x1.0p0;
		else small=0x1.8p0;
		...

here
	i == 10
	x == *d%i == 5 == i/2
but the half-way case does not trigger because z-d == 2 instead of 1
and z[-1] == 0 which should not happen here

Re: printf issues

[Rich Felker]

On Mon, Apr 07, 2014 at 05:36:11PM +0200, Szabolcs Nagy wrote:
> * Morten Welinder <mwelinder@gmail.com> [2014-04-07 10:13:26 -0400]:
> > The frequency of this problem is something like 1 in 5e6.
> > Observations:
> > 
> > *  I only seem to be able to trigger it for %g even though all my
> >    samples print in "e" form.
> > 
> > * The numbers are all roughly the same size: 1e15
> > 
> > * The numbers are all integers ending in 05.  (Except two cases
> >    where the 05 is followed by zeros.)
> > 
> > * The precision is always just below the value that would have
> >   make an exact representation.
> 
> printf("%.12g\n", 1000000000005.0);
> printf("%.11g\n", 500000000045.0);
> printf("%.11g\n", 275000000025.0);
> 
> prints
> 
> 1.00000000001e+12
> 5.0000000005e+11
> 2.7500000003e+11
> 
> in fmt_fp
> 	if (x || d+1!=z) {
> 		long double round = CONCAT(0x1p,LDBL_MANT_DIG);
> 		long double small;
> 		if (*d/i & 1) round += 2;
> 		if (x<i/2) small=0x0.8p0;
> 		else if (x==i/2 && d+1==z) small=0x1.0p0;
> 		else small=0x1.8p0;
> 		...
> 
> here
> 	i == 10
> 	x == *d%i == 5 == i/2
> but the half-way case does not trigger because z-d == 2 instead of 1
> and z[-1] == 0 which should not happen here

Thanks for making the analysis needed to fix this. I've committed the
fix. Hopefully this is the last of such bugs. Perhaps we should devise
a stress test with random inputs and assertions to search for other
bugs. Some ideas that come to mind:

1. Printing both full-precision and roundings to 0...25 places and
   asserting that the roundings are correct based on the full-prec.

2. For full precision outputs, asserting that sum(digits)%3 ==
   numerator%3 for diadic rationals of the form numerator/denominator.

3. Round tripping with strtold.

4. Assertions about formatting such as lack of trailing zeros.

...?

Rich

Re: printf issues

[Szabolcs Nagy]

* Morten Welinder <mwelinder@gmail.com> [2014-04-04 16:01:08 -0400]:
> In fmt_fmt, the rounding decision is done using this test:
> 
>             /* Decide whether to round by probing round+small */
>             if (round+small != round) { ...
> 
> Why is this done with long double?

to support FLT_EVAL_METHOD==2 (for i386 on gcc without sse)

this means that all fp arithmetics is evaluated at long double precision
on the fpu and even if we use different types or cast we have to take
care of excess precision related issues (double-rounding, and various
bugs in older gcc) which may or may not be easy

in this case it seems to me that single precision float can be used
there, we just have to make sure that the excess precision is rounded
away correctly before the round+small != round comparision

> The reason I ask is that the Valgrind situation improves a lot if
> this is done with doubles.
> 
> (Valgrind situation: Valgrind emulates long doubles, poorly, by using
> simple doubles.  See, for example, https://bugs.kde.org/show_bug.cgi?id=164298)
> 

yes, that is known and unfortunate limitation, it's in their manual:
http://valgrind.org/docs/manual/manual-core.html#manual-core.limits

we cannot do much about it: with such broken arithmetics
arbitrarily bad things can happen even if we remove
long doubles from the printf code

Re: printf issues

[Rich Felker]

On Fri, Apr 04, 2014 at 04:01:08PM -0400, Morten Welinder wrote:
> In fmt_fmt, the rounding decision is done using this test:
> 
>             /* Decide whether to round by probing round+small */
>             if (round+small != round) { ...
> 
> Why is this done with long double?
> 
> The reason I ask is that the Valgrind situation improves a lot if
> this is done with doubles.
> 
> (Valgrind situation: Valgrind emulates long doubles, poorly, by using
> simple doubles.  See, for example, https://bugs.kde.org/show_bug.cgi?id=164298)

This is a known issue that needs to be fixed in valgrind. It's
impossible to do anything useful with rounding on x86 with types other
than double, due to excess precision (FLT_EVAL_METHOD==2). This is why
musl uses long double internally everywhere that rounding semantics
matter.

Rich

Re: printf issues

[Morten Welinder]

> [...] excess precision (FLT_EVAL_METHOD==2). This is why
> musl uses long double internally everywhere that rounding semantics
> matter.

That's what I thought, but it's not actually what I see over in src/math/.

If I look in src/math/floor.c I see an explicit cast from double to double
used to get rid of excess precision.  The similar thing ought to work in
fmt_fp.

Morten

On Fri, Apr 4, 2014 at 5:02 PM, Rich Felker <dalias@aerifal.cx> wrote:
> On Fri, Apr 04, 2014 at 04:01:08PM -0400, Morten Welinder wrote:
>> In fmt_fmt, the rounding decision is done using this test:
>>
>>             /* Decide whether to round by probing round+small */
>>             if (round+small != round) { ...
>>
>> Why is this done with long double?
>>
>> The reason I ask is that the Valgrind situation improves a lot if
>> this is done with doubles.
>>
>> (Valgrind situation: Valgrind emulates long doubles, poorly, by using
>> simple doubles.  See, for example, https://bugs.kde.org/show_bug.cgi?id=164298)
>
> This is a known issue that needs to be fixed in valgrind. It's
> impossible to do anything useful with rounding on x86 with types other
> than double, due to excess precision (FLT_EVAL_METHOD==2). This is why
> musl uses long double internally everywhere that rounding semantics
> matter.
>
> Rich

Re: printf issues

[Rich Felker]

On Fri, Apr 04, 2014 at 10:08:47PM -0400, Morten Welinder wrote:
> > [...] excess precision (FLT_EVAL_METHOD==2). This is why
> > musl uses long double internally everywhere that rounding semantics
> > matter.
> 
> That's what I thought, but it's not actually what I see over in src/math/.

I guess I should elaborate that I meant everywhere in the code that I
wrote, which doesn't include anything in src/math except asm.

> If I look in src/math/floor.c I see an explicit cast from double to double
> used to get rid of excess precision.  The similar thing ought to work in
> fmt_fp.

Yes, I think this works, but it's fairly fragile under the possibility
of compiler bugs. FWIW, floor, etc. all have asm versions on i386 so
the excess precision issue doesn't come into play unless you go out of
your way to remove the asm. (This reminds me -- I want to eventually
separate mandatory asm from optimization asm, to make it easier to
test C code that would otherwise be shadowed by asm.)

Rich

Re: printf issues

[Szabolcs Nagy]

* Rich Felker <dalias@aerifal.cx> [2014-04-04 22:50:10 -0400]:
> On Fri, Apr 04, 2014 at 10:08:47PM -0400, Morten Welinder wrote:
> > > [...] excess precision (FLT_EVAL_METHOD==2). This is why
> > > musl uses long double internally everywhere that rounding semantics
> > > matter.
> > 
> > That's what I thought, but it's not actually what I see over in src/math/.
> 
> I guess I should elaborate that I meant everywhere in the code that I
> wrote, which doesn't include anything in src/math except asm.

long double only arith or assuming broken arith would not work for
src/math (it's enough headache to work around the lack of compiler
support for c99 fenv)

the rounding functions (floor etc) are recent code and representative:
- smaller and faster than the original fdlibm/freebsd code
- the only ugliness is fenv related gcc/clang bug workaround
- otherwise clean c99 code relying on c99 + annex f semantics
- not much inline comments (because that would be very repetitive
across the >200 math functions, instead math quirks are documented
on the wiki (for now))

(fwiw they work for i386 as well, the only required asm is sqrt*.s,
removing the rest makes the math code bigger and a bit less precise
but not broken, and yes on old gcc this involves a bit of faith..)


historical notes for the curious:

musl originally used fdlibm with some cleanups, but that was incomplete

so freebsd libm code (and some code from openbsd) was ported to musl
(that was the most well maintained fdlibm version at the time, other
projects did the same: bionic uses it without much modification, julia
lang for scientific computing uses openlibm that is based on the same,
etc.)

unlike freebsd, linux uses the i386 fpu in extended precision mode so
the codebase had to be audited for i386 specific isssues..
(but the freebsd code needed a lot of cleanups and some bug fixes anyway)

the gcc bugs:
- not dropping excess precision at assignment/cast
(infamous gcc bug 323, in freebsd this is worked around by the
STRICT_ASSIGN macro using volatile assignment)
- spuriously dropping excess precision when spilling fpu registers
(annoying because of unpredictable and inconsistent results, but this
did not really come up in the math code because that already stored
most intermediate results and expected that the precision is dropped)
- incorrect rounding of decimal floating-point constants
(not really an issue for math: hexfloats are used where decimal const
would be inexact)
- incorrect handling of long double constants
(freebsd specific, does not matter on linux, but there are ugly
workarounds for it in the freebsd code)

non-gcc issue:
- double-rounding (rounding to long double then to double gives
different results than just one rounding to double)

a clean workaround for all the issues is always using long double vars
(or sometimes float_t/double_t). printf can do this, but in math code
some results need to be correctly rounded to double or float precision
and excess precision is harmful and hard to get rid of.

most of the bugs were fixed in gcc-4.5 (-fexcess-precision=standard or
-std=c99, broken by default) and earlier gccs have -ffloat-store which
fixes the issues to some extent with a lot of extra load/stores.
(clang does not support FLT_EVAL_METHOD!=0, uses sse2 arith on x86)

instead assuming broken excess precision handling and relying on
volatile hacks for correctness it was decided that libm will require
c99 semantics (using the flags above for gcc)

so i removed all the STRICT_ASSIGN macros

however freebsd libm uses float/double instead of float_t/double_t even
if the excess precision is not harmful, so with strict c99 semantics
there were many unwanted load/stores (== slower and less precise code
than necessary on i386, this was fixed, but there are still places that
could be improved)

a few more chapters in the i386 excess precision saga:
- there was no simple and correct sqrt implementation for i386 on
linux, Rich came up with a solution (now it is in glibc as well)
- fma had to be rewritten for i386 (fma and sqrt are the only non-exact
functions that need to be correctly rounded)
- underflow exception may be incorrectly omitted because of double
rounding (i think freebsd still suffers from this), this required
workarounds in asm and c (volatile hacks, many functions were affected)
- gcc -std=c99 drops excess precision at returns (which is incorrect in
c99 but correct in c11) adding useless load/stores
- using sqrtl instead of sqrt/sqrtf internally can be better sometimes
(i did this for acosh where it is provably a strict improvement)
- there is a minor gcc issue about special constants in non-nearest
rounding mode (it recognises pi and uses fldpi instruction)
- posix bug report about M_ macros in math.h when FLT_EVAL_METHOD!=0
- several pending defect reports against iso c about FLT_EVAL_METHOD!=0

> > If I look in src/math/floor.c I see an explicit cast from double to double
> > used to get rid of excess precision.  The similar thing ought to work in
> > fmt_fp.
> 
> Yes, I think this works, but it's fairly fragile under the possibility
> of compiler bugs. FWIW, floor, etc. all have asm versions on i386 so
> the excess precision issue doesn't come into play unless you go out of
> your way to remove the asm. (This reminds me -- I want to eventually
> separate mandatory asm from optimization asm, to make it easier to
> test C code that would otherwise be shadowed by asm.)
> 
> Rich

Re: printf issues

[Rich Felker]

On Fri, Apr 04, 2014 at 08:54:13PM +0200, Szabolcs Nagy wrote:
> > #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
> > long double modfl(long double x, long double *iptr)
> > {
> > return modf(x, (double *)iptr);
> > }
> 
> yes, this is an aliasing violation, nice catch
> 
> the original idea was to allow tail call opt for these wrappers,
> so they are a single branch instruction, we should fix it but
> i think we can rely on that the ptr representations are the same:
> 
> long double modfl(long double x, long double *iptr)
> {
> 	union {long double *ld; double *d;} u = {iptr};
> 	return modf(x, u.d);
> }

No, that doesn't change anything. The problem is that the object being
accessed is aliased with another type, which is illegal. The pointer
type/cast are not a problem, so this "workaround" is working around
the naive outward view of the issue ("bad pointer cast") and not the
actual issue (aliasing).

Rich

Re: printf issues

[Szabolcs Nagy]

* Rich Felker <dalias@aerifal.cx> [2014-04-04 17:00:11 -0400]:
> On Fri, Apr 04, 2014 at 08:54:13PM +0200, Szabolcs Nagy wrote:
> > > #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
> > > long double modfl(long double x, long double *iptr)
> > > {
> > > return modf(x, (double *)iptr);
> > > }
> > 
> > yes, this is an aliasing violation, nice catch
> > 
> > the original idea was to allow tail call opt for these wrappers,
> > so they are a single branch instruction, we should fix it but
> > i think we can rely on that the ptr representations are the same:
> > 
> > long double modfl(long double x, long double *iptr)
> > {
> > 	union {long double *ld; double *d;} u = {iptr};
> > 	return modf(x, u.d);
> > }
> 
> No, that doesn't change anything. The problem is that the object being
> accessed is aliased with another type, which is illegal. The pointer
> type/cast are not a problem, so this "workaround" is working around
> the naive outward view of the issue ("bad pointer cast") and not the
> actual issue (aliasing).

ah yes, my bad

but with temporaries gcc-4.8 does not figure out the tco..
and i don't see other ways around
(other than cheating and adding the single jmp solution to
all arch in asm)

i think only modfl and sincosl are affected so uselessly
shuffling things around the stack is probably not too
expensive

Re: printf issues

[Rich Felker]

On Fri, Apr 04, 2014 at 01:42:30PM -0400, Morten Welinder wrote:
> > I _would_ like this code to be easily adaptable for use outside libc
> > if somebody wants it
> 
> FYI, I have been doing just that for Gnumeric in a variant that always
> rounds ties away from zero.  Two changes would help with making
> the code fit seamlessly into other environments.
> 
> 1. Make "i" in fmt_fp unsigned.  It's used in connection with
>     unsigned values only.

I don't object, but what's the motivation?

> 2. Make "char *s" used to hold "NAN" etc. "const char *s".

Indeed, this is harmless. I assume it's to satisfy that ugly option
that changes the type of string literals from char[] to const char[]?
BTW that option actually just got a lot more problematic with C11; now
that C has _Generic, it potentially changes the semantics of a program
rather than just helping generate warnings.

> Neither of these should make any difference in what the function
> actually does.
> 
> I have run tens of millions random numbers through this function
> looking for differences between it and glibc.  The extra 0s from "%g"
> is the only problem observed.

Nice. I was actually thinking of some numerical tests we could run on
huge random samples, for instance using using theorems about
properties of the digits of the decimal expansion (e.g. n%3==0 iff the
sum of the digits mod 3==0; this works for diadic-rational 'multiples
of 3' too).

> It looks like the LDBL_EPSILON version could be used in
> 
>     roundl.c
>     modfl.c
>     ceill.c
>     floorl.c
> 
> in the definition of TOINT instead of enumerating choices for
> LDBL_MANT_DIG.  It's basically the same thing going on
> there.

Indeed, this would be a welcome change.

> While I was looking for that, I noticed that this modfl fallback looks
> problematic.  Even if long double and double are the same thing
> under the hood, I don't think you can cast pointers like that and
> assume it works.  It needs a temporary.
> 
> #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
> long double modfl(long double x, long double *iptr)
> {
> return modf(x, (double *)iptr);
> }

Agreed. This is UB (an aliasing violation) and should be fixed even if
it makes the function a few bytes larger/a few cycles slower.

Rich