[<prev] [next>] [<thread-prev] [thread-next>] [day] [month] [year] [list]
Message-ID: <20140406230740.GL3034@port70.net>
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.
