This patch removes the large range reduction code and defers to the huge range
reduction code. The first level range reducer supports inputs up to 2^27,
which is way too large given that inputs for sin/cos are typically small
(< 10), and optimizing for a smaller range would give a significant speedup.
Input values above 2^27 are practically never used, so there is no reason for
supporting range reduction between 2^27 and 2^48. Removing it significantly
simplifies code and enables further speedups. There is about a 2.3x slowdown
in this range due to __branred being extremely slow (a better algorithm could
easily more than double performance).
* sysdeps/ieee754/dbl-64/s_sin.c (reduce_sincos_2): Remove function.
(do_sincos_2): Likewise.
(__sin): Remove middle range reduction case.
(__cos): Likewise.
* sysdeps/ieee754/dbl-64/s_sincos.c (__sincos): Remove middle range
reduction case.
This series of patches removes the slow patchs from sin, cos and sincos.
Besides greatly simplifying the implementation, the new version is also much
faster for inputs up to PI (41% faster) and for large inputs needing range
reduction (27% faster).
ULP is ~0.55 with no errors found after testing 1.6 billion inputs across most
of the range with mpsin and mpcos. The number of incorrectly rounded results
(ie. ULP >0.5) is at most ~2750 per million inputs between 0.125 and 0.5,
the average is ~850 per million between 0 and PI.
Tested on AArch64 and x86_64 with no regressions.
The first patch removes the slow paths for the cases where the input is small
and doesn't require range reduction. Update ULP tables for sin, cos and sincos
on AArch64 and x86_64.
* sysdeps/aarch64/libm-test-ulps: Update ULP for sin, cos, sincos.
* sysdeps/ieee754/dbl-64/s_sin.c (__sin): Remove slow paths for small
inputs.
(__cos): Likewise.
* sysdeps/x86_64/fpu/libm-test-ulps: Update ULP for sin, cos, sincos.
Continuing the move of libm aliases to common macros that can create
_FloatN / _FloatNx aliases in future, this patch converts some dbl-64
functions to using libm_alias_double, thereby eliminating the need for
some ldbl-opt wrappers.
This patch deliberately limits what functions are converted so that it
can be verified by comparison of stipped binaries. Specifically, atan
and tan are excluded because they first need converting to being weak
aliases; fma is omitted as it has additional complications with
versions in other directories (removing the ldbl-opt version can
e.g. cause the ldbl-128 version to be used instead of dbl-64); and
functions that have both dbl-64/wordsize-64 and ldbl-opt versions are
excluded because ldbl-opt currently always wraps dbl-64 function
versions, so changing those will result in platforms using both
ldbl-opt and dbl-64/wordsize-64 (i.e. alpha) starting to use the
dbl-64/wordsize-64 versions of those functions (which is good, as an
optimization, but still best separated from the present patch to get
better validation).
Tested for x86_64, and tested with build-many-glibcs.py that installed
stripped shared libraries are unchanged by the patch.
* sysdeps/ieee754/dbl-64/s_asinh.c: Include <libm-alias-double.h>.
(asinh): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_cbrt.c: Include <libm-alias-double.h>.
(cbrt): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_copysign.c: Include
<libm-alias-double.h>.
(copysign): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_erf.c: Include <libm-alias-double.h>.
(erf): Define using libm_alias_double.
(erfc): Likewise.
* sysdeps/ieee754/dbl-64/s_expm1.c: Include <libm-alias-double.h>.
(expm1): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_fabs.c: Include <libm-alias-double.h>.
(fabs): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_fromfp.c (fromfp): Define using
libm_alias_double.
* sysdeps/ieee754/dbl-64/s_fromfp_main.c: Include
<libm-alias-double.h>.
* sysdeps/ieee754/dbl-64/s_fromfpx.c (fromfpx): Define using
libm_alias_double.
* sysdeps/ieee754/dbl-64/s_getpayload.c: Include
<libm-alias-double.h>.
(getpayload): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_llrint.c: Include
<libm-alias-double.h>.
(llrint): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_lrint.c: Include <libm-alias-double.h>.
(lrint): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_nextup.c: Include
<libm-alias-double.h>.
(nextup): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_roundeven.c: Include
<libm-alias-double.h>.
(roundeven): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_setpayload.c (setpayload): Define using
libm_alias_double.
* sysdeps/ieee754/dbl-64/s_setpayload_main.c: Include
<libm-alias-double.h>.
* sysdeps/ieee754/dbl-64/s_setpayloadsig.c (setpayloadsig): Define
using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_sin.c: Include <libm-alias-double.h>.
(cos): Define using libm_alias_double.
(sin): Likewise.
* sysdeps/ieee754/dbl-64/s_sincos.c: Include
<libm-alias-double.h>.
(sincos): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_tanh.c: Include <libm-alias-double.h>.
(tanh): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_totalorder.c: Include
<libm-alias-double.h>.
(totalorder): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_totalordermag.c: Include
<libm-alias-double.h>.
(totalordermag): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/s_ufromfp.c (ufromfp): Define using
libm_alias_double.
* sysdeps/ieee754/dbl-64/s_ufromfpx.c (ufromfpx): Define using
libm_alias_double.
* sysdeps/ieee754/dbl-64/wordsize-64/s_getpayload.c: Include
<libm-alias-double.h>.
(getpayload): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/wordsize-64/s_roundeven.c: Include
<libm-alias-double.h>.
(roundeven): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/wordsize-64/s_setpayload_main.c: Include
<libm-alias-double.h>.
* sysdeps/ieee754/dbl-64/wordsize-64/s_totalorder.c: Include
<libm-alias-double.h>.
(totalorder): Define using libm_alias_double.
* sysdeps/ieee754/dbl-64/wordsize-64/s_totalordermag.c: Include
<libm-alias-double.h>.
(totalordermag): Define using libm_alias_double.
* sysdeps/ieee754/ldbl-opt/s_copysign.c (copysignl): Only define
libc compat symbol here.
* sysdeps/ieee754/ldbl-opt/s_asinh.c: Remove file.
* sysdeps/ieee754/ldbl-opt/s_cbrt.c: Likewise.
* sysdeps/ieee754/ldbl-opt/s_erf.c: Likewise.
* sysdeps/ieee754/ldbl-opt/s_expm1.c: Likewise.
* sysdeps/ieee754/ldbl-opt/s_fabs.c: Likewise.
* sysdeps/ieee754/ldbl-opt/s_llrint.c: Likewise.
* sysdeps/ieee754/ldbl-opt/s_lrint.c: Likewise.
* sysdeps/ieee754/ldbl-opt/s_sin.c: Likewise.
* sysdeps/ieee754/ldbl-opt/s_sincos.c: Likewise.
* sysdeps/ieee754/ldbl-opt/s_tanh.c: Likewise.
Like the previous change, make the quadrant shift a boolean to make it
clearer that we will do at most a single rotation of the quadrants to
compute the cosine from the sine function.
This does not affect codegen.
For k1 in 1 and 3, n can only have values of 0 and 2, so checking k1 &
2 is equivalent to checking n & 2. We prefer the latter so that we
don't use k1 for anything other than selecting the quadrant in
do_sincos_1, thus dropping it completely.
The previous logic was:
"Compute sine for the value and based on the new rotated quadrant
(k1) negate the value if we're in the fourth quadrant."
With this change, the logic now is:
"Compute sine for the value and negate it if we were either (1) in
the fourth quadrant or (2) we actually wanted the cosine and were
in the third quadrant."
* sysdeps/ieee754/dbl-64/s_sin.c (do_sincos_1): Check N
instead of K1.
The do_sincos_* functions are helpers to compute sin/cos, where they
get cosine by computing sine for the next quadrant. This is decided
with the value of K passed to it, which is the amount by which to
shift the quadrant. Since we will only need the shift to be 0 or 1,
we make K a bool to make that explicit.
* sysdeps/ieee754/dbl-64/s_sin.c (do_sincos_1): Rename K to
SHIFT_QUADRANT and make it bool.
(do_sincos_2): Likewise.
(sloww): Likewise.
(sloww1): Likewise.
(__sin): Adjust calls to do_sincos_1 and do_sincos_2.
(__cos): Likewise.
These are remaining cases where we can deduce and conclude that the
sign of the result should be the same as the sign of the input being
checked. For example, for sin(x), the sign of the result is the same
as the result itself for x < pi. Likewise, for sine values where x
after range reduction falls into this range and its sign is preserved.
* sysdeps/ieee754/dbl-64/s_sin.c (do_sincos_1): Use copysign
instead of ternary condition.
(do_sincos_2): Likewise.
(__sin): Likewise.
(__cos): Likewise.
(slow): Likewise.
(sloww): Likewise.
(sloww1): Likewise.
(bsloww): Likewise.
(bsloww1): Likewise.
This is the first very simple substitution of ternary conditions for
correction adjustments with __copysign for positive constants.
* sysdeps/ieee754/dbl-64/s_sin.c (do_cos_slow): use copysign
instead of ternary condition.
(do_sin_slow): Likewise.
(do_sincos_1): Likewise.
(do_sincos_2): Likewise.
(__cos): Likewise.
(sloww): Likewise.
(sloww1): Likewise.
(sloww2): Likewise.
(bsloww): Likewise.
(bsloww1): Likewise.
(bsloww2): Likewise.
Simplify the code a bit by consolidating sign checks in slow1 and
slow2 into __sin at the higher level.
* sysdeps/ieee754/dbl-64/s_sin.c (slow1): Consolidate sign
check from here...
(slow2): ... and here...
(__sin): ... to here.
The support functions for sin and cos have a lot of identical
functionality, so inlining them gives a pretty decent jump in
functionality: ~19% in the sincos function. On SPEC2006 this
translates to about 2.1% in the tonto test.
* sysdeps/ieee754/dbl-64/s_sin.c (do_cos): Mark as inline.
(do_cos_slow): Likewise.
(do_sin): Likewise.
(do_sin_slow): Likewise.
(slow): Likewise.
(slow1): Likewise.
(slow2): Likewise.
(sloww): Likewise.
(sloww1): Likewise.
(sloww2): Likewise.
(bsloww): Likewise.
(bsloww1): Likewise.
(bsloww2): Likewise.
(cslow2): Likewise.
The only code looks slightly different from do_sin but on closer
examination, should give exactly the same result. Drop it in favour
of the do_sin function call.
* sysdeps/ieee754/dbl-64/s_sin.c (__sin): Use do_sin.
All calls to do_cos are preceded by code that partitions x into a
larger double that gives an offset into the sincos table and a smaller
double that is used in a polynomial computation. Consolidate all of
them into do_cos and do_sin to reduce code duplication.
* sysdeps/ieee754/dbl-64/s_sin.c (do_cos): Accept X and DX as input
arguments. Consolidate input partitioning from callers here.
(do_cos_slow): Likewise.
(do_sin): Likewise.
(do_sin_slow): Likewise.
(do_sincos_1): Remove the no longer necessary input partitioning.
(do_sincos_2): Likewise.
(__sin): Likewise.
(__cos): Likewise.
(slow1): Likewise.
(slow2): Likewise.
(sloww1): Likewise.
(sloww2): Likewise.
(bsloww1): Likewise.
(bsloww2): Likewise.
(cslow2): Likewise.
The sin and cos code is inconsistent about its use of fabs to get the
absolute value of X where in some places it conditionalizes the code
while in others it uses fabs. fabs seems to be a better candidate in
most cases because it avoids a branch. Similarly there is an attempt
to make it easier for the compiler to emit conditional assignment
instructions (like fcsel on aarch64) where it can, by isolating
conditional assignment constructs from the rest of the expression.
A further benefit of this change is to identify common constructs
across functions and consolidate them in future patches.
* sysdeps/ieee754/dbl-64/s_sin.c (do_cos_slow): Use ternary
instead of if/else.
(do_sin_slow): Likewise.
(do_sincos_1): Use fabs instead of if/else.
(do_sincos_2): Likewise.
(__sin): Likewise.
(__cos): Likewise.
(slow2): Likewise.
(sloww): Likewise.
(sloww1): Likewise. Drop argument M.
(sloww2): Use fabs instead of if/else.
(bsloww): Likewise.
(bsloww1): Likewise.
(bsloww2): Likewise.
This patch reshuffles the reduce_and_compute code so that the
structure matches other code structures of the same type elsewhere in
s_sin.c and s_sincos.c. This is the beginning of an attempt to
consolidate and reduce code duplication in functions in s_sin.c to
make it easier to read and possibly also easier for the compiler to
optimize.
* sysdeps/ieee754/dbl-64/s_sin.c (reduce_and_compute):
Consolidate switch cases 0 and 2.
During the sincos consolidation I made two mistakes, one was a logical
error due to which cos(0x1.8475e5afd4481p+0) returned
sin(0x1.8475e5afd4481p+0) instead.
The second issue was an error in negating inputs for the correct
quadrants for sine. I could not find a suitable test case for this
despite running a program to search for such an input for a couple of
hours.
Following patch fixes both issues. Tested on x86_64. Thanks to Matt
Clay for identifying the issue.
[BZ #20357]
* sysdeps/ieee754/dbl-64/s_sin.c (sloww): Fix up condition
to call __mpsin/__mpcos and to negate values.
* math/auto-libm-test-in: Add test.
* math/auto-libm-test-out: Regenerate.
Like the previous change, exploit the fact that computation for sin
and cos is identical except that it is apart by a quadrant. Also
remove csloww, csloww1 and csloww2 since they can easily be expressed
in terms of sloww, sloww1 and sloww2.
The sin and cos computation for this range of input is identical
except for a difference in quadrants by 1. Exploit that fact and the
common argument reduction to reduce computations for sincos.
Range reduction needs to be done only once for sin and cos, so copy
over all of the relevant functions (__sin, __cos, reduce_and_compute)
and consolidate common code.
Include the __sin and __cos functions as local static copies to allow
deper optimization of the functions. This change shows an improvement
of about 17% in the min case and 12.5% in the mean case for the sincos
microbenchmark on x86_64.
* sysdeps/ieee754/dbl-64/s_sin.c (__sin)[IN_SINCOS]: Mark function
static and don't set or restore rounding.
(__cos)[IN_SINCOS]: Likewise.
* sysdeps/ieee754/dbl-64/s_sincos.c: Include s_sin.c.
(__sincos): Set and restore rounding mode. Remove check for infinite
or NaN input.
Similar to various other bugs in this area, some sin and sincos
implementations do not raise the underflow exception for subnormal
arguments, when the result is tiny and inexact. This patch forces the
exception in a similar way to previous fixes.
Tested for x86_64, x86, mips64 and powerpc.
[BZ #16526]
[BZ #16538]
* sysdeps/ieee754/dbl-64/s_sin.c: Include <float.h>.
(__sin): Force underflow exception for arguments with small
absolute value.
* sysdeps/ieee754/flt-32/k_sinf.c: Include <float.h>.
(__kernel_sinf): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128/k_sincosl.c: Include <float.h>.
(__kernel_sincosl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128/k_sinl.c: Include <float.h>.
(__kernel_sinl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128ibm/k_sincosl.c: Include <float.h>.
(__kernel_sincosl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-128ibm/k_sinl.c: Include <float.h>.
(__kernel_sinl): Force underflow exception for arguments with
small absolute value.
* sysdeps/ieee754/ldbl-96/k_sinl.c: Include <float.h>.
(__kernel_sinl): Force underflow exception for arguments with
small absolute value.
* sysdeps/powerpc/fpu/k_sinf.c: Include <float.h>.
(__kernel_sinf): Force underflow exception for arguments with
small absolute value.
* math/auto-libm-test-in: Add more tests of sin and sincos.
* math/auto-libm-test-out: Regenerated.
In 84ba214c, I removed some redundant sign computations and in the
process, I incorrectly got rid of a temporary variable, thus passing
the absolute value of the input to bsloww1. This caused #16623.
This fix undoes the incorrect change.
This patch consolidates the multiple copies of code that looks up sin
and cos of a number from the lookup table and computes the final
value, into static functions. This does not have a noticeable
performance impact since the functions are inlined by gcc.
There is further scope for consolidation in the functions but they
cause a more noticable impact on performance (>5%) due to which I have
held back on them.
Removed more redundant computations in the slow paths of the sin and
cos functions. The notable change is the passing of the most
significant bits of X to the slow functions to check if X is positive
so that just the absolute value of x can be passed and the repeated
ABS() operation is avoided.
There are multiple points in the code where the absolute value of a
number is computed multiple times or is computed even though the value
can only be positive. This change removes those redundant
computations. Tested on x86_64 to verify that there were no
regressions in the testsuite.
- Remove redundant mynumber union definitions
- Clean up a clumsy ternary operator
- Rename TAYLOR_SINCOS to TAYLOR_SIN since we're only expanding the
sin Taylor series in it.
