Date: Mon, 7 Apr 2014 01:07:41 +0200 From: Szabolcs Nagy <nsz@...t70.net> To: musl@...ts.openwall.com Subject: Re: printf issues * Rich Felker <dalias@...ifal.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
Powered by blists - more mailing lists
Confused about mailing lists and their use? Read about mailing lists on Wikipedia and check out these guidelines on proper formatting of your messages.