* printf issues @ 2014-04-04 13:32 Morten Welinder 2014-04-04 14:12 ` Rich Felker 2014-04-04 14:15 ` Szabolcs Nagy 0 siblings, 2 replies; 27+ messages in thread From: Morten Welinder @ 2014-04-04 13:32 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 13:32 printf issues Morten Welinder @ 2014-04-04 14:12 ` Rich Felker 2014-04-04 14:15 ` Szabolcs Nagy 1 sibling, 0 replies; 27+ messages in thread From: Rich Felker @ 2014-04-04 14:12 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 13:32 printf issues Morten Welinder 2014-04-04 14:12 ` Rich Felker @ 2014-04-04 14:15 ` Szabolcs Nagy 2014-04-04 14:35 ` Morten Welinder 1 sibling, 1 reply; 27+ messages in thread From: Szabolcs Nagy @ 2014-04-04 14:15 UTC (permalink / raw) To: musl * 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) ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 14:15 ` Szabolcs Nagy @ 2014-04-04 14:35 ` Morten Welinder 2014-04-04 14:56 ` Szabolcs Nagy 2014-04-04 15:07 ` Rich Felker 0 siblings, 2 replies; 27+ messages in thread From: Morten Welinder @ 2014-04-04 14:35 UTC (permalink / raw) To: musl > 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) ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 14:35 ` Morten Welinder @ 2014-04-04 14:56 ` Szabolcs Nagy 2014-04-04 15:07 ` Rich Felker 1 sibling, 0 replies; 27+ messages in thread From: Szabolcs Nagy @ 2014-04-04 14:56 UTC (permalink / raw) To: musl * 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 14:35 ` Morten Welinder 2014-04-04 14:56 ` Szabolcs Nagy @ 2014-04-04 15:07 ` Rich Felker 2014-04-04 17:42 ` Morten Welinder 1 sibling, 1 reply; 27+ messages in thread From: Rich Felker @ 2014-04-04 15:07 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 15:07 ` Rich Felker @ 2014-04-04 17:42 ` Morten Welinder 2014-04-04 18:54 ` Szabolcs Nagy 2014-04-04 20:58 ` Rich Felker 0 siblings, 2 replies; 27+ messages in thread From: Morten Welinder @ 2014-04-04 17:42 UTC (permalink / raw) To: musl > 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 17:42 ` Morten Welinder @ 2014-04-04 18:54 ` Szabolcs Nagy 2014-04-04 20:01 ` Morten Welinder 2014-04-04 21:00 ` Rich Felker 2014-04-04 20:58 ` Rich Felker 1 sibling, 2 replies; 27+ messages in thread From: Szabolcs Nagy @ 2014-04-04 18:54 UTC (permalink / raw) To: musl * 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); } ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 18:54 ` Szabolcs Nagy @ 2014-04-04 20:01 ` Morten Welinder 2014-04-04 20:22 ` Morten Welinder ` (2 more replies) 2014-04-04 21:00 ` Rich Felker 1 sibling, 3 replies; 27+ messages in thread From: Morten Welinder @ 2014-04-04 20:01 UTC (permalink / raw) To: musl 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); > } ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 20:01 ` Morten Welinder @ 2014-04-04 20:22 ` Morten Welinder 2014-04-04 21:08 ` Rich Felker 2014-04-04 20:54 ` Szabolcs Nagy 2014-04-04 21:02 ` Rich Felker 2 siblings, 1 reply; 27+ messages in thread From: Morten Welinder @ 2014-04-04 20:22 UTC (permalink / raw) To: musl 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); >> } ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 20:22 ` Morten Welinder @ 2014-04-04 21:08 ` Rich Felker 2014-04-04 22:50 ` Morten Welinder 0 siblings, 1 reply; 27+ messages in thread From: Rich Felker @ 2014-04-04 21:08 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 21:08 ` Rich Felker @ 2014-04-04 22:50 ` Morten Welinder 2014-04-05 0:01 ` Morten Welinder 0 siblings, 1 reply; 27+ messages in thread From: Morten Welinder @ 2014-04-04 22:50 UTC (permalink / raw) To: musl > 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 22:50 ` Morten Welinder @ 2014-04-05 0:01 ` Morten Welinder 2014-04-05 1:41 ` Rich Felker 2014-04-07 7:29 ` Rich Felker 0 siblings, 2 replies; 27+ messages in thread From: Morten Welinder @ 2014-04-05 0:01 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-05 0:01 ` Morten Welinder @ 2014-04-05 1:41 ` Rich Felker 2014-04-07 7:29 ` Rich Felker 1 sibling, 0 replies; 27+ messages in thread From: Rich Felker @ 2014-04-05 1:41 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-05 0:01 ` Morten Welinder 2014-04-05 1:41 ` Rich Felker @ 2014-04-07 7:29 ` Rich Felker 2014-04-07 13:40 ` Morten Welinder 1 sibling, 1 reply; 27+ messages in thread From: Rich Felker @ 2014-04-07 7:29 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-07 7:29 ` Rich Felker @ 2014-04-07 13:40 ` Morten Welinder 2014-04-07 14:13 ` Morten Welinder 0 siblings, 1 reply; 27+ messages in thread From: Morten Welinder @ 2014-04-07 13:40 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-07 13:40 ` Morten Welinder @ 2014-04-07 14:13 ` Morten Welinder 2014-04-07 15:36 ` Szabolcs Nagy 0 siblings, 1 reply; 27+ messages in thread From: Morten Welinder @ 2014-04-07 14:13 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-07 14:13 ` Morten Welinder @ 2014-04-07 15:36 ` Szabolcs Nagy 2014-04-07 18:04 ` Rich Felker 0 siblings, 1 reply; 27+ messages in thread From: Szabolcs Nagy @ 2014-04-07 15:36 UTC (permalink / raw) To: musl * 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-07 15:36 ` Szabolcs Nagy @ 2014-04-07 18:04 ` Rich Felker 0 siblings, 0 replies; 27+ messages in thread From: Rich Felker @ 2014-04-07 18:04 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 20:01 ` Morten Welinder 2014-04-04 20:22 ` Morten Welinder @ 2014-04-04 20:54 ` Szabolcs Nagy 2014-04-04 21:02 ` Rich Felker 2 siblings, 0 replies; 27+ messages in thread From: Szabolcs Nagy @ 2014-04-04 20:54 UTC (permalink / raw) To: musl * 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 20:01 ` Morten Welinder 2014-04-04 20:22 ` Morten Welinder 2014-04-04 20:54 ` Szabolcs Nagy @ 2014-04-04 21:02 ` Rich Felker 2014-04-05 2:08 ` Morten Welinder 2 siblings, 1 reply; 27+ messages in thread From: Rich Felker @ 2014-04-04 21:02 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 21:02 ` Rich Felker @ 2014-04-05 2:08 ` Morten Welinder 2014-04-05 2:50 ` Rich Felker 0 siblings, 1 reply; 27+ messages in thread From: Morten Welinder @ 2014-04-05 2:08 UTC (permalink / raw) To: musl > [...] 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-05 2:08 ` Morten Welinder @ 2014-04-05 2:50 ` Rich Felker 2014-04-06 23:07 ` Szabolcs Nagy 0 siblings, 1 reply; 27+ messages in thread From: Rich Felker @ 2014-04-05 2:50 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-05 2:50 ` Rich Felker @ 2014-04-06 23:07 ` Szabolcs Nagy 0 siblings, 0 replies; 27+ messages in thread From: Szabolcs Nagy @ 2014-04-06 23:07 UTC (permalink / raw) To: musl * 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 18:54 ` Szabolcs Nagy 2014-04-04 20:01 ` Morten Welinder @ 2014-04-04 21:00 ` Rich Felker 2014-04-04 21:10 ` Szabolcs Nagy 1 sibling, 1 reply; 27+ messages in thread From: Rich Felker @ 2014-04-04 21:00 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 21:00 ` Rich Felker @ 2014-04-04 21:10 ` Szabolcs Nagy 0 siblings, 0 replies; 27+ messages in thread From: Szabolcs Nagy @ 2014-04-04 21:10 UTC (permalink / raw) To: musl * 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
* Re: printf issues 2014-04-04 17:42 ` Morten Welinder 2014-04-04 18:54 ` Szabolcs Nagy @ 2014-04-04 20:58 ` Rich Felker 1 sibling, 0 replies; 27+ messages in thread From: Rich Felker @ 2014-04-04 20:58 UTC (permalink / raw) To: musl 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 ^ permalink raw reply [flat|nested] 27+ messages in thread
end of thread, other threads:[~2014-04-07 18:04 UTC | newest] Thread overview: 27+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 2014-04-04 13:32 printf issues Morten Welinder 2014-04-04 14:12 ` Rich Felker 2014-04-04 14:15 ` Szabolcs Nagy 2014-04-04 14:35 ` Morten Welinder 2014-04-04 14:56 ` Szabolcs Nagy 2014-04-04 15:07 ` Rich Felker 2014-04-04 17:42 ` Morten Welinder 2014-04-04 18:54 ` Szabolcs Nagy 2014-04-04 20:01 ` Morten Welinder 2014-04-04 20:22 ` Morten Welinder 2014-04-04 21:08 ` Rich Felker 2014-04-04 22:50 ` Morten Welinder 2014-04-05 0:01 ` Morten Welinder 2014-04-05 1:41 ` Rich Felker 2014-04-07 7:29 ` Rich Felker 2014-04-07 13:40 ` Morten Welinder 2014-04-07 14:13 ` Morten Welinder 2014-04-07 15:36 ` Szabolcs Nagy 2014-04-07 18:04 ` Rich Felker 2014-04-04 20:54 ` Szabolcs Nagy 2014-04-04 21:02 ` Rich Felker 2014-04-05 2:08 ` Morten Welinder 2014-04-05 2:50 ` Rich Felker 2014-04-06 23:07 ` Szabolcs Nagy 2014-04-04 21:00 ` Rich Felker 2014-04-04 21:10 ` Szabolcs Nagy 2014-04-04 20:58 ` Rich Felker
Code repositories for project(s) associated with this public inbox https://git.vuxu.org/mirror/musl/ This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox; as well as URLs for NNTP newsgroup(s).