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, 79 insertions, 0 deletions

diff --git a/math/sincosf_data.c b/math/sincosf_data.c
new file mode 100644
index 0000000..ee6aeb9
--- /dev/null
+++ b/math/sincosf_data.c
@@ -0,0 +1,79 @@
+/*
+ * Data definition for sinf, cosf and sincosf.
+ *
+ * Copyright (c) 2018, Arm Limited.
+ * SPDX-License-Identifier: Apache-2.0
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#if !WANT_SINGLEPREC
+
+#include <stdint.h>
+#include <math.h>
+#include "math_config.h"
+#include "sincosf.h"
+
+/* The constants and polynomials for sine and cosine.  The 2nd entry
+   computes -cos (x) rather than cos (x) to get negation for free.  */
+sincos_t sincosf_table[2] =
+{
+  {
+    { 1.0, -1.0, -1.0, 1.0 },
+#if HAVE_FAST_ROUND
+    0x1.45F306DC9C883p-1,
+#else
+    0x1.45F306DC9C883p+23,
+#endif
+    0x1.921FB54442D18p0,
+    0x1p0,
+    -0x1.ffffffd0c621cp-2,
+    0x1.55553e1068f19p-5,
+    -0x1.6c087e89a359dp-10,
+    0x1.99343027bf8c3p-16,
+    -0x1.555545995a603p-3,
+    0x1.1107605230bc4p-7,
+    -0x1.994eb3774cf24p-13
+  },
+  {
+    { 1.0, -1.0, -1.0, 1.0 },
+#if HAVE_FAST_ROUND
+    0x1.45F306DC9C883p-1,
+#else
+    0x1.45F306DC9C883p+23,
+#endif
+    0x1.921FB54442D18p0,
+    -0x1p0,
+    0x1.ffffffd0c621cp-2,
+    -0x1.55553e1068f19p-5,
+    0x1.6c087e89a359dp-10,
+    -0x1.99343027bf8c3p-16,
+    -0x1.555545995a603p-3,
+    0x1.1107605230bc4p-7,
+    -0x1.994eb3774cf24p-13
+  }
+};
+
+/* Table with 4/PI to 192 bit precision.  To avoid unaligned accesses
+   only 8 new bits are added per entry, making the table 4 times larger.  */
+const uint32_t inv_pio4[24] =
+{
+  0xa2,       0xa2f9,	  0xa2f983,   0xa2f9836e,
+  0xf9836e4e, 0x836e4e44, 0x6e4e4415, 0x4e441529,
+  0x441529fc, 0x1529fc27, 0x29fc2757, 0xfc2757d1,
+  0x2757d1f5, 0x57d1f534, 0xd1f534dd, 0xf534ddc0,
+  0x34ddc0db, 0xddc0db62, 0xc0db6295, 0xdb629599,
+  0x6295993c, 0x95993c43, 0x993c4390, 0x3c439041
+};
+
+#endif