Improve performance of sinf/cosf/sincosf

This patch is a complete rewrite of sinf, cosf and sincosf.  The new version
is significantly faster, as well as simple and accurate.
The worst-case ULP is 0.56072, maximum relative error is 0.5303p-23 over all
4 billion inputs.  In non-nearest rounding modes the error is 1ULP.

The algorithm uses 3 main cases: small inputs which don't need argument
reduction, small inputs which need a simple range reduction and large inputs
requiring complex range reduction.  The code uses approximate integer
comparisons to quickly decide between these cases - on some targets this may
be slow, so this can be configured to use floating point comparisons.

The small range reducer uses a single reduction step to handle values up to
120.0.  It is fastest on targets which support inlined round instructions.

The large range reducer uses integer arithmetic for simplicity.  It does a
32x96 bit multiply to compute a 64-bit modulo result.  This is more than
accurate enough to handle the worst-case cancellation for values close to
an integer multiple of PI/4.  It could be further optimized, however it is
already much faster than necessary.

1 files changed, 3 insertions, 0 deletions

diff --git a/math/include/mathlib.h b/math/include/mathlib.h
index 43962ab..dff4922 100644
--- a/math/include/mathlib.h
+++ b/math/include/mathlib.h
@@ -25,3 +25,6 @@ float exp2f(float);
 float logf(float);
 float log2f(float);
 float powf(float, float);
+float sinf(float);
+float cosf(float);
+void sincosf(float, float*, float*);