the volatile hack in STRICT_ASSIGN is only needed if
assignment is not respected and excess precision is kept.
gcc -fexcess-precision=standard and -ffloat-store both
respect assignment and musl use these flags by default.
i kept the macro for now so the workaround may be used
for bad compilers in the future.
old code was correct only if the result was stored (without the
excess precision) or musl was compiled with -ffloat-store.
now we use STRICT_ASSIGN to work around the issue.
(see note 160 in c11 section 6.8.6.4)
old code was correct only if the result was stored (without the
excess precision) or musl was compiled with -ffloat-store.
(see note 160 in n1570.pdf section 6.8.6.4)
The long double adjustment was wrong:
The usual check is
mant_bits & 0x7ff == 0x400
before doing a mant_bits++ or mant_bits-- adjustment since
this is the only case when rounding an inexact ld80 into
double can go wrong. (only in nearest rounding mode)
After such a check the ++ and -- is ok (the mantissa will end
in 0x401 or 0x3ff).
fma is a bit different (we need to add 3 numbers with correct
rounding: hi_xy + lo_xy + z so we should survive two roundings
at different places without precision loss)
The adjustment in fma only checks for zero low bits
mant_bits & 0x3ff == 0
this way the adjusted value is correct when rounded to
double or *less* precision.
(this is an important piece in the fma puzzle)
Unfortunately in this case the -- is not a correct adjustment
because mant_bits might underflow so further checks are needed
and this was the source of the bug.
old: 2*atan2(sqrt(1-x),sqrt(1+x))
new: atan2(fabs(sqrt((1-x)*(1+x))),x)
improvements:
* all edge cases are fixed (sign of zero in downward rounding)
* a bit faster (here a single call is about 131ns vs 162ns)
* a bit more precise (at most 1ulp error on 1M uniform random
samples in [0,1), the old formula gave some 2ulp errors as well)
the old formula atan2(1,sqrt((1+x)/(1-x))) was faster but
could give nan result at x=1 when the rounding mode is
FE_DOWNWARD (so 1-1 == -0 and 2/-0 == -inf), the new formula
gives -0 at x=+-1 with downward rounding.
old code saved/restored the fenv (the new code is only as slow
as that when inexact is not set before the call, but some other
flag is set and the rounding is inexact, which is rare)
before:
bench_nearbyint_exact 5000000 N 261 ns/op
bench_nearbyint_inexact_set 5000000 N 262 ns/op
bench_nearbyint_inexact_unset 5000000 N 261 ns/op
after:
bench_nearbyint_exact 10000000 N 94.99 ns/op
bench_nearbyint_inexact_set 25000000 N 65.81 ns/op
bench_nearbyint_inexact_unset 10000000 N 94.97 ns/op
zero, one, two, half are replaced by const literals
The policy was to use the f suffix for float consts (1.0f),
but don't use suffix for long double consts (these consts
can be exactly represented as double).
Underflow exception is only raised when the result is
invalid, but fmod is always exact. x87 has a denormalization
exception, but that's nonstandard. And the superflous *1.0
will be optimized away by any compiler that does not honor
signaling nans.
Some code assumed ldexp(x, 1) is faster than 2.0*x,
but ldexp is a wrapper around scalbn which uses
multiplications inside, so this optimization is
wrong.
This commit also fixes fmal which accidentally
used ldexp instead of ldexpl loosing precision.
There are various additional changes from the
work-in-progress const cleanups.
Some long double consts were stored in two doubles as a workaround
for x86_64 and i386 with the following comment:
/* Long double constants are slow on these arches, and broken on i386. */
This is most likely old gcc bug related to the default x87 fpu
precision setting (it's double instead of double extended on BSD).
A faster workaround for spurious inexact exceptions
when the result cannot be represented. The old code
actually could be wrong, because gcc reordered the
integer conversion and the exception check.
Note that the new fesetround has slightly different semantics:
Storing the floating-point environment with fnstenv makes the
next fldenv (or fldcw) "non-signaling", so unmasked and pending
exceptions does not invoke the exception handler.
(These are rare since exceptions are handled immediately and by
default all exceptions are masked anyway. But if one manually
unmasks an exception in the control word then either sets the
corresponding exception flag in the status word or the execution
of an exception raising floating-point operation gets interrupted
then it may happen).
So the old implementation did not trap in some rare cases
where the new implementation traps.
However POSIX does not specify anything like the x87 exception
handling traps and the fnstenv/fldenv pair is significantly slower
than the fnstcw/fldcw pair (new code is about 5x faster here and
it's dominated by the function call overhead).
Szabolcs Nagy