- fma(x, y, z) returns z, not NaN, if z is infinite, x and y are finite,
x*y overflows, and x*y and z have opposite signs.
- fma(x, y, z) doesn't generate an overflow, underflow, or inexact exception
if z is NaN or infinite, as per IEEE 754R.
- If the rounding mode is set to FE_DOWNWARD, fma(1.0, 0.0, -0.0) is -0.0,
not +0.0.
This makes little difference in float precision, but in double
precision gives a speedup of about 30% on amd64 (A64 CPU) and i386
(A64). This depends on fabs[f]() being inline and efficient. The
bit fiddling (or any use of SET_HIGH_WORD(), which libm does too
much because it was best on old 32-bit machines) always causes
packing overheads and sometimes causes stalls in the packing, since
it operates on only part of a variable in the double precision case.
It apparently did cause stalls in a critical path here.
fabs(x+0.0)+fabs(y+0.0) when mixing NaNs. This improves
consistency of the result by making it harder for the compiler to reorder
the operands. (FP addition is not necessarily commutative because the
order of operands makes a difference on some machines iff the operands are
both NaNs.)
e_rem_pio2.c:
This case goes up to about 2**20pi/2, but the comment about it said that
it goes up to about 2**19pi/2.
It went too far above 2**pi/2, giving a multiplier fn with 21 significant
bits in some cases. This would be harmful except for a numerical
accident. It happens that the terms of the approximation to pi/2,
when rounded to 33 bits so that multiplications by 20-bit fn's are
exact, happen to be rounded to 32 bits so multiplications by 21-bit
fn's are exact too, so the bug only complicates the error analysis (we
might lose a bit of accuracy but have bits to spare).
e_rem_pio2f.c:
The bogus comment in e_rem_pio2.c was copied and the code was changed
to be bug-for-bug compatible with it, except the limit was made 90
ulps smaller than necessary. The approximation to pi/2 was not
modified except for discarding some of it.
The same rough error analysis that justifies the limit of 2**20pi/2
for double precision only justifies a limit of 2**18pi/2 for float
precision. We depended on exhaustive testing to check the magic numbers
for float precision. More exaustive testing shows that we can go up
to 2**28pi/2 using a 53+25 bit approximation to pi/2 for float precision,
with a the maximum error for cosf() and sinf() unchanged at 0.5009
ulps despite the maximum error in rem_pio2f being ~0.25 ulps. Implement
this.
gives an average speedup of about 12 cycles or 17% for
9pi/4 < |x| <= 2**19pi/2 and a smaller speedup for larger x, and a
small speeddown for |x| <= 9pi/4 (only 1-2 cycles average, but that
is 4%).
Inlining this is less likely to bust caches than inlining the float
version since it is much smaller (about 220 bytes text and rodata) and
has many fewer branches. However, the float version was already large
due to its manual inlining of the branches and also the polynomial
evaluations.
__kernel_rem_pio2(). This simplifies analysis of aliasing and thus
results in better code for the usual case where __kernel_rem_pio2()
is not called. In particular, when __ieee854_rem_pio2[f]() is inlined,
it normally results in y[] being returned in registers. I couldn't
get this to work using the restrict qualifier.
In float precision, this saves 2-3% in most cases on amd64 and i386
(A64) despite it not being inlined in float precision yet. In double
precision, this has high variance, with an average gain of 2% for
amd64 and 0.7% for i386 (but a much larger gain for usual cases) and
some losses.
this function and its callers cosf(), sinf() and tanf() don't waste time
converting values from doubles to floats and back for |x| > 9pi/4.
All these functions were optimized a few years ago to mostly use doubles
internally and across the __kernel*() interfaces but not across the
__ieee754_rem_pio2f() interface.
This saves about 40 cycles in cosf(), sinf() and tanf() for |x| > 9pi/4
on amd64 (A64), and about 20 cycles on i386 (A64) (except for cosf()
and sinf() in the upper range). 40 cycles is about 35% for |x| < 9pi/4
<= 2**19pi/2 and about 5% for |x| > 2**19pi/2. The saving is much
larger on amd64 than on i386 since the conversions are not easy to
optimize except on i386 where some of them are automatic and others
are optimized invalidly. amd64 is still about 10% slower in cosf()
and tanf() in the lower range due to conversion overhead.
This also gives a tiny speedup for |x| <= 9pi/4 on amd64 (by simplifying
the code). It also avoids compiler bugs and/or additional slowness
in the conversions on (not yet supported) machines where double_t !=
double.
e_rem_pio2.c:
Float and double precision didn't work because init_jk[] was 1 too small.
It needs to be 2 larger than you might expect, and 1 larger than it was
for these precisions, since its test for recomputing needs a margin of
47 bits (almost 2 24-bit units).
init_jk[] seems to be barely enough for extended and quad precisions.
This hasn't been completely verified. Callers now get about 24 bits
of extra precision for float, and about 19 for double, but only about
8 for extended and quad. 8 is not enough for callers that want to
produce extra-precision results, but current callers have rounding
errors of at least 0.8 ulps, so another 1/2**8 ulps of error from the
reduction won't affect them much.
Add a comment about some of the magic for init_jk[].
e_rem_pio2.c:
Double precision worked in practice because of a compensating off-by-1
error here. Extended precision was asked for, and it executed exactly
the same code as the unbroken double precision.
e_rem_pio2f.c:
Float precision worked in practice because of a compensating off-by-1
error here. Double precision was asked for, and was almost needed,
since the cosf() and sinf() callers want to produce extra-precision
results, at least internally so that their error is only 0.5009 ulps.
However, the extra precision provided by unbroken float precision is
enough, and the double-precision code has extra overheads, so the
off-by-1 error cost about 5% in efficiency on amd64 and i386.
variations (e500 currently), this provides a gcc-level FPU emulation and is an
alternative approach to the recently introduced kernel-level emulation
(FPU_EMU).
Approved by: cognet (mentor)
MFp4: e500
fabs(), a conditional branch, and sign adjustments of 3 variables for
x < 0 when the branch is taken. In double precision, even when the
branch is perfectly predicted, this saves about 10 cycles or 10% on
amd64 (A64) and i386 (A64) for the negative half of the range, but
makes little difference for the positive half of the range. In float
precision, it also saves about 4 cycles for the positive half of the
range on i386, and many more cycles in both halves on amd64 (28 in the
negative half and 11 in the positive half for tanf), but the amd64
times for float precision are anomalously slow so the larger
improvement is only a side effect.
Previous commits arranged for the x < 0 case to be handled simply:
- one part of the rounding method uses the magic number 0x1.8p52
instead of the usual 0x1.0p52. The latter is required for large |x|,
but it doesn't work for negative x and we don't need it for large |x|.
- another part of the rounding method no longer needs to add `half'.
It would have needed to add -half for negative x.
- removing the "quick check no cancellation" in the double precision
case removed the need to take the absolute value of the quadrant
number.
Add my noncopyright in e_rem_pio2.c
FP-to-FP method to round to an integer on all arches, and convert this
to an int using FP-to-integer conversion iff irint() is not available.
This is cleaner and works well on at least ia64, where it saves 20-30
cycles or about 10% on average for 9Pi/4 < |x| <= 32pi/2 (should be
similar up to 2**19pi/2, but I only tested the smaller range).
After the previous commit to e_rem_pio2.c removed the "quick check no
cancellation" non-optimization, the result of the FP-to-integer
conversion is not needed so early, so using irint() became a much
smaller optimization than when it was committed.
An earlier commit message said that cos, cosf, sin and sinf were equally
fast on amd64 and i386 except for cos and sin on i386. Actually, cos
and sin on amd64 are equally fast to cosf and sinf on i386 (~88 cycles),
while cosf and sinf on amd64 are not quite equally slow to cos and sin
on i386 (average 115 cycles with more variance).
9pi/2 < |x| < 32pi/2 since it is only a small or negative optimation
and it gets in the way of further optimizations. It did one more
branch to avoid some integer operations and to use a different
dependency on previous results. The branches are fairly predictable
so they are usually not a problem, so whether this is a good
optimization depends mainly on the timing for the previous results,
which is very machine-dependent. On amd64 (A64), this "optimization"
is a pessimization of about 1 cycle or 1%; on ia64, it is an
optimization of about 2 cycles or 1%; on i386 (A64), it is an
optimization of about 5 cycles or 4%; on i386 (Celeron P2) it is an
optimization of about 4 cycles or 3% for cos but a pessimization of
about 5 cycles for sin and 1 cycle for tan. I think the new i386
(A64) slowness is due to an pipeline stall due to an avoidable
load-store mismatch (so the old timing was better), and the i386
(Celeron) variance is due to its branch predictor not being too good.
the the double to int conversion operation which is very slow on these
arches. Assume that the current rounding mode is the default of
round-to-nearest and use rounding operations in this mode instead of
faking this mode using the round-towards-zero mode for conversion to
int. Round the double to an integer as a double first and as an int
second since the double result is needed much earler.
Double rounding isn't a problem since we only need a rough approximation.
We didn't support other current rounding modes and produce much larger
errors than before if called in a non-default mode.
This saves an average about 10 cycles on amd64 (A64) and about 25 on
i386 (A64) for x in the above range. In some cases the saving is over
25%. Most cases with |x| < 1000pi now take about 88 cycles for cos
and sin (with certain CFLAGS, etc.), except on i386 where cos and sin
(but not cosf and sinf) are much slower at 111 and 121 cycles respectivly
due to the compiler only optimizing well for float precision. A64
hardware cos and sin are slower at 105 cycles on i386 and 110 cycles
on amd64.
the same as lrint() except it returns int instead of long. Though the
extern lrint() is fairly fast on these arches, it still takes about
12 cycles longer than the inline version, and 12 cycles is a lot in
applications where [li]rint() is used to avoid slow conversions that
are only a couple of times slower.
This is only for internal use. The libm versions of *rint*() should
also be inline, but that would take would take more header engineering.
Implementing irint() instead of lrint() also avoids a conflict with
the extern declaration of the latter.
on i386 (A64), 5 cycles on amd64 (A64), and 3 cycles on ia64). gcc
tends to generate very bad code for accessing floating point values
as bits except when the integer accesses have the same width as the
floating point values, and direct accesses to bit-fields (as is common
only for long double precision) always gives such accesses. Use the
expsign access method, which is good for 80-bit long doubles and
hopefully no worse for 128-bit long doubles. Now the generated code
is less bad. There is still unnecessary copying of the arg on amd64
and i386 and mysterious extra slowness on amd64.
pi/4 <= |x| <= 3pi/4. Use the same branch ladder as for float precision.
Remove the optimization for |x| near pi/2 and don't do it near the
multiples of pi/2 in the newly optimized range, since it requires
fairly large code to handle only relativley few cases. Ifdef out
optimization for |x| <= pi/4 since this case can't occur because it
is done in callers.
On amd64 (A64), for cos() and sin() with uniformly distributed args,
no cache misses, some parallelism in the caller, and good but not great
CC and CFLAGS, etc., this saves about 40 cycles or 38% in the newly
optimized range, or about 27% on average across the range |x| <= 2pi
(~65 cycles for most args, while the A64 hardware fcos and fsin take
~75 cycles for half the args and 125 cycles for the other half). The
speedup for tan() is much smaller, especially relatively. The speedup
on i386 (A64) is slightly smaller, especially relatively. i386 is
still much slower than amd64 here (unlike in the float case where it
is slightly faster).
saves an average of about 8 cycles or 5% on A64 (amd64 and i386 --
more in cycles but about the same percentage on i386, and more with
old versions of gcc) with good CFLAGS and some parallelism in the
caller. As usual, it takes a couple more multiplications so it will
be slower on old machines.
Convert to __FBSDID().
optimization of about 10% for cos(x), sin(x) and tan(x) on
|x| < 2**19*pi/2. We didn't do this before because __ieee754__rem_pio2()
is too large and complicated for gcc-3.3 to inline very well. We don't
do this for float precision because it interferes with optimization
of the usual (?) case (|x| < 9pi/4) which is manually inlined for float
precision only.
This has some rough edges:
- some static data is duplicated unnecessarily. There isn't much after
the recent move of large tables to k_rem_pio2.c, and some static data
is duplicated to good affect (all the data static const, so that the
compiler can evaluate expressions like 2*pio2 at compile time and
generate even more static data for the constant for this).
- extern inline is used (for the same reason as in previous inlining of
k_cosf.c etc.), but C99 apparently doesn't allow extern inline
functions with static data, and gcc will eventually warn about this.
Convert to __FBSDID().
Indent __ieee754_rem_pio2()'s declaration consistently (its style was
made inconsistent with fdlibm a while ago, so complete this).
Fix __ieee754_rem_pio2()'s return type to match its prototype. Someone
changed too many ints to int32_t's when fixing the assumption that all
ints are int32_t's.
