This is the mail archive of the
newlib@sourceware.org
mailing list for the newlib project.
Re: [PATCH] Improve performance of sinf/cosf/sincosf
Here is the v3 version:
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.
Simple benchmark showing speedup factor on AArch64 for various ranges:
range 0.7853982 sinf 1.7 cosf 2.2 sincosf 2.8
range 1.570796 sinf 1.9 cosf 1.9 sincosf 2.7
range 3.141593 sinf 2.0 cosf 2.0 sincosf 3.5
range 6.283185 sinf 2.3 cosf 2.3 sincosf 4.2
range 125.6637 sinf 2.9 cosf 3.0 sincosf 5.1
range 1.1259e15 sinf 26.8 cosf 26.8 sincosf 45.2
ChangeLog:
2018-05-18 Wilco Dijkstra <wdijkstr@arm.com>
* newlib/libm/common/Makefile.in: Regenerated.
* newlib/libm/common/Makefile.am: Add sinf.c, cosf.c, sincosf.c
sincosf.h, sincosf_data.c. Add -fbuiltin -fno-math-errno to CFLAGS.
* newlib/libm/common/math_config.h: Add HAVE_FAST_ROUND, HAVE_FAST_LROUND,
roundtoint, converttoint, force_eval_float, force_eval_double, eval_as_float,
eval_as_double, likely, unlikely.
* newlib/libm/common/cosf.c: New file.
* newlib/libm/common/sinf.c: Likewise.
* newlib/libm/common/sincosf.h: Likewise.
* newlib/libm/common/sincosf.c: Likewise.
* newlib/libm/common/sincosf_data.c: Likewise.
* newlib/libm/math/sf_cos.c: Add #if to build conditionally.
* newlib/libm/math/sf_sin.c: Likewise.
* newlib/libm/math/wf_sincos.c: Likewise.
--
diff --git a/newlib/libm/common/Makefile.am b/newlib/libm/common/Makefile.am
index 3e2620c7ebd27fe09b74695cd377a983aa5cc584..423d82eb79152c11040ecfe42338f85b3f73f38b 100644
--- a/newlib/libm/common/Makefile.am
+++ b/newlib/libm/common/Makefile.am
@@ -25,7 +25,7 @@ fsrc = sf_finite.c sf_copysign.c sf_modf.c sf_scalbn.c \
sf_scalbln.c sf_trunc.c \
sf_exp.c sf_exp2.c sf_exp2_data.c sf_log.c sf_log_data.c \
sf_log2.c sf_log2_data.c sf_pow_log2_data.c sf_pow.c \
- math_errf.c
+ sf_sin.c sf_cos.c sf_sincos.c sf_sincos_data.c math_errf.c
lsrc = atanl.c cosl.c sinl.c tanl.c tanhl.c frexpl.c modfl.c ceill.c fabsl.c \
floorl.c log1pl.c expm1l.c acosl.c asinl.c atan2l.c coshl.c sinhl.c \
@@ -38,6 +38,7 @@ lsrc = atanl.c cosl.c sinl.c tanl.c tanhl.c frexpl.c modfl.c ceill.c fabsl.c \
sl_finite.c
libcommon_la_LDFLAGS = -Xcompiler -nostdlib
+lib_a_CFLAGS = -fbuiltin -fno-math-errno
if USE_LIBTOOL
noinst_LTLIBRARIES = libcommon.la
@@ -52,7 +53,7 @@ lib_a_SOURCES = $(src) $(fsrc)
if HAVE_LONG_DOUBLE
lib_a_SOURCES += $(lsrc)
endif # HAVE_LONG_DOUBLE
-lib_a_CFLAGS = $(AM_CFLAGS)
+lib_a_CFLAGS += $(AM_CFLAGS)
noinst_DATA =
endif # USE_LIBTOOL
diff --git a/newlib/libm/common/Makefile.in b/newlib/libm/common/Makefile.in
index 6a0ba73d24d039cc610d417df474a51ccba223b9..e4fb7c62f09e25388e61e8acd7296c6197cb399c 100644
--- a/newlib/libm/common/Makefile.in
+++ b/newlib/libm/common/Makefile.in
@@ -55,6 +55,7 @@ build_triplet = @build@
host_triplet = @host@
@HAVE_LONG_DOUBLE_TRUE@@USE_LIBTOOL_TRUE@am__append_1 = $(lsrc)
@HAVE_LONG_DOUBLE_TRUE@@USE_LIBTOOL_FALSE@am__append_2 = $(lsrc)
+@USE_LIBTOOL_FALSE@am__append_3 = $(AM_CFLAGS)
DIST_COMMON = $(srcdir)/../../Makefile.shared $(srcdir)/Makefile.in \
$(srcdir)/Makefile.am
subdir = common
@@ -114,6 +115,8 @@ am__objects_2 = lib_a-sf_finite.$(OBJEXT) lib_a-sf_copysign.$(OBJEXT) \
lib_a-sf_log.$(OBJEXT) lib_a-sf_log_data.$(OBJEXT) \
lib_a-sf_log2.$(OBJEXT) lib_a-sf_log2_data.$(OBJEXT) \
lib_a-sf_pow_log2_data.$(OBJEXT) lib_a-sf_pow.$(OBJEXT) \
+ lib_a-sf_sin.$(OBJEXT) lib_a-sf_cos.$(OBJEXT) \
+ lib_a-sf_sincos.$(OBJEXT) lib_a-sf_sincos_data.$(OBJEXT) \
lib_a-math_errf.$(OBJEXT)
am__objects_3 = lib_a-atanl.$(OBJEXT) lib_a-cosl.$(OBJEXT) \
lib_a-sinl.$(OBJEXT) lib_a-tanl.$(OBJEXT) \
@@ -169,7 +172,8 @@ am__objects_6 = sf_finite.lo sf_copysign.lo sf_modf.lo sf_scalbn.lo \
sf_nearbyint.lo sf_remquo.lo sf_round.lo sf_scalbln.lo \
sf_trunc.lo sf_exp.lo sf_exp2.lo sf_exp2_data.lo sf_log.lo \
sf_log_data.lo sf_log2.lo sf_log2_data.lo sf_pow_log2_data.lo \
- sf_pow.lo math_errf.lo
+ sf_pow.lo sf_sin.lo sf_cos.lo sf_sincos.lo sf_sincos_data.lo \
+ math_errf.lo
am__objects_7 = atanl.lo cosl.lo sinl.lo tanl.lo tanhl.lo frexpl.lo \
modfl.lo ceill.lo fabsl.lo floorl.lo log1pl.lo expm1l.lo \
acosl.lo asinl.lo atan2l.lo coshl.lo sinhl.lo expl.lo \
@@ -361,7 +365,7 @@ fsrc = sf_finite.c sf_copysign.c sf_modf.c sf_scalbn.c \
sf_scalbln.c sf_trunc.c \
sf_exp.c sf_exp2.c sf_exp2_data.c sf_log.c sf_log_data.c \
sf_log2.c sf_log2_data.c sf_pow_log2_data.c sf_pow.c \
- math_errf.c
+ sf_sin.c sf_cos.c sf_sincos.c sf_sincos_data.c math_errf.c
lsrc = atanl.c cosl.c sinl.c tanl.c tanhl.c frexpl.c modfl.c ceill.c fabsl.c \
floorl.c log1pl.c expm1l.c acosl.c asinl.c atan2l.c coshl.c sinhl.c \
@@ -374,6 +378,7 @@ lsrc = atanl.c cosl.c sinl.c tanl.c tanhl.c frexpl.c modfl.c ceill.c fabsl.c \
sl_finite.c
libcommon_la_LDFLAGS = -Xcompiler -nostdlib
+lib_a_CFLAGS = -fbuiltin -fno-math-errno $(am__append_3)
@USE_LIBTOOL_TRUE@noinst_LTLIBRARIES = libcommon.la
@USE_LIBTOOL_TRUE@libcommon_la_SOURCES = $(src) $(fsrc) \
@USE_LIBTOOL_TRUE@ $(am__append_1)
@@ -381,7 +386,6 @@ libcommon_la_LDFLAGS = -Xcompiler -nostdlib
@USE_LIBTOOL_TRUE@noinst_DATA = objectlist.awk.in
@USE_LIBTOOL_FALSE@noinst_LIBRARIES = lib.a
@USE_LIBTOOL_FALSE@lib_a_SOURCES = $(src) $(fsrc) $(am__append_2)
-@USE_LIBTOOL_FALSE@lib_a_CFLAGS = $(AM_CFLAGS)
#
# documentation rules
@@ -944,6 +948,30 @@ lib_a-sf_pow.o: sf_pow.c
lib_a-sf_pow.obj: sf_pow.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_pow.obj `if test -f 'sf_pow.c'; then $(CYGPATH_W) 'sf_pow.c'; else $(CYGPATH_W) '$(srcdir)/sf_pow.c'; fi`
+lib_a-sf_sin.o: sf_sin.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_sin.o `test -f 'sf_sin.c' || echo '$(srcdir)/'`sf_sin.c
+
+lib_a-sf_sin.obj: sf_sin.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_sin.obj `if test -f 'sf_sin.c'; then $(CYGPATH_W) 'sf_sin.c'; else $(CYGPATH_W) '$(srcdir)/sf_sin.c'; fi`
+
+lib_a-sf_cos.o: sf_cos.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_cos.o `test -f 'sf_cos.c' || echo '$(srcdir)/'`sf_cos.c
+
+lib_a-sf_cos.obj: sf_cos.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_cos.obj `if test -f 'sf_cos.c'; then $(CYGPATH_W) 'sf_cos.c'; else $(CYGPATH_W) '$(srcdir)/sf_cos.c'; fi`
+
+lib_a-sf_sincos.o: sf_sincos.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_sincos.o `test -f 'sf_sincos.c' || echo '$(srcdir)/'`sf_sincos.c
+
+lib_a-sf_sincos.obj: sf_sincos.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_sincos.obj `if test -f 'sf_sincos.c'; then $(CYGPATH_W) 'sf_sincos.c'; else $(CYGPATH_W) '$(srcdir)/sf_sincos.c'; fi`
+
+lib_a-sf_sincos_data.o: sf_sincos_data.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_sincos_data.o `test -f 'sf_sincos_data.c' || echo '$(srcdir)/'`sf_sincos_data.c
+
+lib_a-sf_sincos_data.obj: sf_sincos_data.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sf_sincos_data.obj `if test -f 'sf_sincos_data.c'; then $(CYGPATH_W) 'sf_sincos_data.c'; else $(CYGPATH_W) '$(srcdir)/sf_sincos_data.c'; fi`
+
lib_a-math_errf.o: math_errf.c
$(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-math_errf.o `test -f 'math_errf.c' || echo '$(srcdir)/'`math_errf.c
diff --git a/newlib/libm/common/math_config.h b/newlib/libm/common/math_config.h
index c5667844a85d025f0d0173e82faca54d24a3d632..6d87376ab5c3ac88fe0982977660f1e27b416dde 100644
--- a/newlib/libm/common/math_config.h
+++ b/newlib/libm/common/math_config.h
@@ -1,5 +1,5 @@
/* Configuration for math routines.
- Copyright (c) 2017 ARM Ltd. All rights reserved.
+ Copyright (c) 2017-2018 Arm Ltd. All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
@@ -43,6 +43,44 @@
# define WANT_ERRNO_UFLOW (WANT_ROUNDING && WANT_ERRNO)
#endif
+/* Compiler can inline round as a single instruction. */
+#ifndef HAVE_FAST_ROUND
+# if __aarch64__
+# define HAVE_FAST_ROUND 1
+# else
+# define HAVE_FAST_ROUND 0
+# endif
+#endif
+
+/* Compiler can inline lround, but not (long)round(x). */
+#ifndef HAVE_FAST_LROUND
+# if __aarch64__ && (100*__GNUC__ + __GNUC_MINOR__) >= 408 && __NO_MATH_ERRNO__
+# define HAVE_FAST_LROUND 1
+# else
+# define HAVE_FAST_LROUND 0
+# endif
+#endif
+
+#if HAVE_FAST_ROUND
+# define TOINT_INTRINSICS 1
+
+static inline double_t
+roundtoint (double_t x)
+{
+ return round (x);
+}
+
+static inline uint64_t
+converttoint (double_t x)
+{
+# if HAVE_FAST_LROUND
+ return lround (x);
+# else
+ return (long) round (x);
+# endif
+}
+#endif
+
#ifndef TOINT_INTRINSICS
# define TOINT_INTRINSICS 0
#endif
@@ -109,12 +147,56 @@ issignalingf_inline (float x)
return 2 * (ix ^ 0x00400000) > 2u * 0x7fc00000;
}
+/* Force the evaluation of a floating-point expression for its side-effect. */
+#if __aarch64__ && __GNUC__
+static inline void
+force_eval_float (float x)
+{
+ __asm__ __volatile__ ("" : "+w" (x));
+}
+static inline void
+force_eval_double (double x)
+{
+ __asm__ __volatile__ ("" : "+w" (x));
+}
+#else
+static inline void
+force_eval_float (float x)
+{
+ volatile float y = x;
+}
+static inline void
+force_eval_double (double x)
+{
+ volatile double y = x;
+}
+#endif
+
+/* Evaluate an expression as the specified type, normally a type
+ cast should be enough, but compilers implement non-standard
+ excess-precision handling, so when FLT_EVAL_METHOD != 0 then
+ these functions may need to be customized. */
+static inline float
+eval_as_float (float x)
+{
+ return x;
+}
+static inline double
+eval_as_double (double x)
+{
+ return x;
+}
+
#ifdef __GNUC__
# define HIDDEN __attribute__ ((__visibility__ ("hidden")))
# define NOINLINE __attribute__ ((noinline))
+# define likely(x) __builtin_expect (!!(x), 1)
+# define unlikely(x) __builtin_expect (x, 0)
#else
# define HIDDEN
# define NOINLINE
+# define likely(x) (x)
+# define unlikely(x) (x)
#endif
HIDDEN float __math_oflowf (unsigned long);
diff --git a/newlib/libm/common/sf_cos.c b/newlib/libm/common/sf_cos.c
new file mode 100644
index 0000000000000000000000000000000000000000..4a59d7d72d4c80d75b97aff7d377101cb456e590
--- /dev/null
+++ b/newlib/libm/common/sf_cos.c
@@ -0,0 +1,88 @@
+/* Single-precision cos function.
+ Copyright (c) 2018 Arm Ltd. All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+ 1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ 2. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+ 3. The name of the company may not be used to endorse or promote
+ products derived from this software without specific prior written
+ permission.
+
+ THIS SOFTWARE IS PROVIDED BY ARM LTD ``AS IS AND ANY EXPRESS OR IMPLIED
+ WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+ MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ IN NO EVENT SHALL ARM LTD BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
+
+#include "fdlibm.h"
+#if !__OBSOLETE_MATH
+
+#include <stdint.h>
+#include <math.h>
+#include "math_config.h"
+#include "sf_sincos.h"
+
+/* Fast cosf implementation. Worst-case ULP is 0.56072, maximum relative
+ error is 0.5303p-23. A single-step signed range reduction is used for
+ small values. Large inputs have their range reduced using fast integer
+ arithmetic.
+*/
+float
+cosf (float y)
+{
+ double x = y;
+ double s;
+ int n;
+ sincos_t *p = &sincosf_table[0];
+
+ if (abstop12 (y) < abstop12 (pio4))
+ {
+ double x2 = x * x;
+
+ if (unlikely (abstop12 (y) < abstop12 (0x1p-12f)))
+ return 1.0f;
+
+ return sinf_poly (x, x2, p, 1);
+ }
+ else if (likely (abstop12 (y) < abstop12 (120.0f)))
+ {
+ x = reduce_fast (x, p, &n);
+
+ /* Setup the signs for sin and cos. */
+ s = p->sign[n & 3];
+
+ if (n & 2)
+ p = &sincosf_table[1];
+
+ return sinf_poly (x * s, x * x, p, n ^ 1);
+ }
+ else if (abstop12 (y) < abstop12 (INFINITY))
+ {
+ uint32_t xi = asuint (y);
+ int sign = xi >> 31;
+
+ x = reduce_large (xi, &n);
+
+ /* Setup signs for sin and cos - include original sign. */
+ s = p->sign[(n + sign) & 3];
+
+ if ((n + sign) & 2)
+ p = &sincosf_table[1];
+
+ return sinf_poly (x * s, x * x, p, n ^ 1);
+ }
+ else
+ return __math_invalidf (y);
+}
+
+#endif
diff --git a/newlib/libm/common/sf_sin.c b/newlib/libm/common/sf_sin.c
new file mode 100644
index 0000000000000000000000000000000000000000..8250a1b4aee0b5db5fb88f947fcd25dde97de777
--- /dev/null
+++ b/newlib/libm/common/sf_sin.c
@@ -0,0 +1,92 @@
+/* Single-precision sin function.
+ Copyright (c) 2018 Arm Ltd. All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+ 1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ 2. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+ 3. The name of the company may not be used to endorse or promote
+ products derived from this software without specific prior written
+ permission.
+
+ THIS SOFTWARE IS PROVIDED BY ARM LTD ``AS IS AND ANY EXPRESS OR IMPLIED
+ WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+ MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ IN NO EVENT SHALL ARM LTD BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
+
+#include "fdlibm.h"
+#if !__OBSOLETE_MATH
+
+#include <math.h>
+#include "math_config.h"
+#include "sf_sincos.h"
+
+/* Fast sinf implementation. Worst-case ULP is 0.56072, maximum relative
+ error is 0.5303p-23. A single-step signed range reduction is used for
+ small values. Large inputs have their range reduced using fast integer
+ arithmetic.
+*/
+float
+sinf (float y)
+{
+ double x = y;
+ double s;
+ int n;
+ sincos_t *p = &sincosf_table[0];
+
+ if (abstop12 (y) < abstop12 (pio4))
+ {
+ s = x * x;
+
+ if (unlikely (abstop12 (y) < abstop12 (0x1p-12f)))
+ {
+ if (unlikely (abstop12 (y) < abstop12 (0x1p-126f)))
+ /* Force underflow for tiny y. */
+ force_eval_float (s);
+ return y;
+ }
+
+ return sinf_poly (x, s, p, 0);
+ }
+ else if (likely (abstop12 (y) < abstop12 (120.0f)))
+ {
+ x = reduce_fast (x, p, &n);
+
+ /* Setup the signs for sin and cos. */
+ s = p->sign[n & 3];
+
+ if (n & 2)
+ p = &sincosf_table[1];
+
+ return sinf_poly (x * s, x * x, p, n);
+ }
+ else if (abstop12 (y) < abstop12 (INFINITY))
+ {
+ uint32_t xi = asuint (y);
+ int sign = xi >> 31;
+
+ x = reduce_large (xi, &n);
+
+ /* Setup signs for sin and cos - include original sign. */
+ s = p->sign[(n + sign) & 3];
+
+ if ((n + sign) & 2)
+ p = &sincosf_table[1];
+
+ return sinf_poly (x * s, x * x, p, n);
+ }
+ else
+ return __math_invalidf (y);
+}
+
+#endif
diff --git a/newlib/libm/common/sf_sincos.h b/newlib/libm/common/sf_sincos.h
new file mode 100644
index 0000000000000000000000000000000000000000..2f3b2f83b6f66dbe351cac0d2d896d928081ff48
--- /dev/null
+++ b/newlib/libm/common/sf_sincos.h
@@ -0,0 +1,172 @@
+/* Header for single-precision sin/cos/sincos functions.
+ Copyright (c) 2018 Arm Ltd. All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+ 1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ 2. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+ 3. The name of the company may not be used to endorse or promote
+ products derived from this software without specific prior written
+ permission.
+
+ THIS SOFTWARE IS PROVIDED BY ARM LTD ``AS IS AND ANY EXPRESS OR IMPLIED
+ WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+ MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ IN NO EVENT SHALL ARM LTD BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
+
+#include <stdint.h>
+#include <math.h>
+#include "math_config.h"
+
+/* PI * 2^-64. */
+static const double pi64 = 0x1.921FB54442D18p-62;
+/* PI / 4. */
+static const double pio4 = 0x1.921FB54442D18p-1;
+
+typedef const struct
+{
+ double sign[4];
+ double hpi_inv, hpi, c0, c1, c2, c3, c4, s1, s2, s3;
+} sincos_t;
+
+extern sincos_t sincosf_table[2] HIDDEN;
+
+extern const uint32_t inv_pio4[] HIDDEN;
+
+/* abstop12 assumes floating point reinterpret is fast by default.
+ If floating point comparisons are faster, define PREFER_FLOAT_COMPARISON. */
+#if PREFER_FLOAT_COMPARISON
+static inline float
+abstop12 (float x)
+{
+ return fabsf (x);
+}
+#else
+static inline uint32_t
+abstop12 (float x)
+{
+ return (asuint (x) >> 20) & 0x7ff;
+}
+#endif
+
+/* Compute the sine and cosine of inputs X and X2 (X squared), using the
+ polynomial P and store the results in SINP and COSP. N is the quadrant,
+ if odd the cosine and sine polynomials are swapped. */
+static inline void
+sincosf_poly (double x, double x2, sincos_t *p, int n, float *sinp, float *cosp)
+{
+ double x3, x4, x5, x6, s, c, c1, c2, s1;
+
+ x4 = x2 * x2;
+ x3 = x2 * x;
+ c2 = p->c3 + x2 * p->c4;
+ s1 = p->s2 + x2 * p->s3;
+
+ /* Swap sin/cos result based on quadrant. */
+ float *tmp = (n & 1 ? cosp : sinp);
+ cosp = (n & 1 ? sinp : cosp);
+ sinp = tmp;
+
+ c1 = p->c0 + x2 * p->c1;
+ x5 = x3 * x2;
+ x6 = x4 * x2;
+
+ s = x + x3 * p->s1;
+ c = c1 + x4 * p->c2;
+
+ *sinp = s + x5 * s1;
+ *cosp = c + x6 * c2;
+}
+
+/* Return the sine of inputs X and X2 (X squared) using the polynomial P.
+ N is the quadrant, and if odd the cosine polynomial is used. */
+static inline float
+sinf_poly (double x, double x2, sincos_t *p, int n)
+{
+ double x3, x4, x6, x7, s, c, c1, c2, s1;
+
+ if ((n & 1) == 0)
+ {
+ x3 = x * x2;
+ s1 = p->s2 + x2 * p->s3;
+
+ x7 = x3 * x2;
+ s = x + x3 * p->s1;
+
+ return s + x7 * s1;
+ }
+ else
+ {
+ x4 = x2 * x2;
+ c2 = p->c3 + x2 * p->c4;
+ c1 = p->c0 + x2 * p->c1;
+
+ x6 = x4 * x2;
+ c = c1 + x4 * p->c2;
+
+ return c + x6 * c2;
+ }
+}
+
+/* Fast range reduction using single multiply-subtract. Return the modulo of
+ X as a value between -PI/4 and PI/4 and store the quadrant in NP.
+ The values for PI/2 and 2/PI are accessed via P. Since PI/2 as a double
+ is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4,
+ only 2 multiplies are required and the result is accurate for |X| <= 120.0.
+ Use round/lround if inlined, otherwise convert to int. To avoid inaccuracies
+ introduced by truncating negative values, compute the quadrant * 2^24. */
+static inline double
+reduce_fast (double x, sincos_t *p, int *np)
+{
+ double r;
+#if TOINT_INTRINSICS
+ r = x * p->hpi_inv;
+ *np = converttoint (r);
+ return x - roundtoint (r) * p->hpi;
+#else
+ r = x * p->hpi_inv;
+ int n = ((int32_t)r + 0x800000) >> 24;
+ *np = n;
+ return x - n * p->hpi;
+#endif
+}
+
+/* Reduce the range of XI to a multiple of PI/4 using fast integer arithmetic.
+ XI is a reinterpreted float and must be >= 2.0f (the sign bit is ignored).
+ Return the modulo between -PI/4 and PI/4 and store the quadrant in NP.
+ Reduction uses a table of 4/PI with 192 bits of precision. A 32x96->128 bit
+ multiply computes the exact 2.62-bit fixed-point modulo. Since the result
+ can have at most 29 leading zeros after the binary point, the double
+ precision result is accurate to 33 bits. */
+static inline double
+reduce_large (uint32_t xi, int *np)
+{
+ const uint32_t *arr = &inv_pio4[(xi >> 26) & 15];
+ int shift = (xi >> 23) & 7;
+ uint64_t n, res0, res1, res2;
+
+ xi = (xi & 0xffffff) | 0x800000;
+ xi <<= shift;
+
+ res0 = xi * arr[0];
+ res1 = (uint64_t)xi * arr[4];
+ res2 = (uint64_t)xi * arr[8];
+ res0 = (res2 >> 32) | (res0 << 32);
+ res0 += res1;
+
+ n = (res0 + (1ULL << 61)) >> 62;
+ res0 -= n << 62;
+ double x = (int64_t)res0;
+ *np = n;
+ return x * pi64;
+}
diff --git a/newlib/libm/common/sf_sincos.c b/newlib/libm/common/sf_sincos.c
new file mode 100644
index 0000000000000000000000000000000000000000..74c9fc4cbcfc0dd1e788266a4b2e1ac7570616ad
--- /dev/null
+++ b/newlib/libm/common/sf_sincos.c
@@ -0,0 +1,104 @@
+/* Single-precision sincos function.
+ Copyright (c) 2018 Arm Ltd. All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+ 1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ 2. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+ 3. The name of the company may not be used to endorse or promote
+ products derived from this software without specific prior written
+ permission.
+
+ THIS SOFTWARE IS PROVIDED BY ARM LTD ``AS IS AND ANY EXPRESS OR IMPLIED
+ WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+ MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ IN NO EVENT SHALL ARM LTD BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
+
+#include "fdlibm.h"
+#if !__OBSOLETE_MATH
+
+#include <stdint.h>
+#include <math.h>
+#include "math_config.h"
+#include "sf_sincos.h"
+
+/* Fast sincosf implementation. Worst-case ULP is 0.56072, maximum relative
+ error is 0.5303p-23. A single-step signed range reduction is used for
+ small values. Large inputs have their range reduced using fast integer
+ arithmetic.
+*/
+void
+sincosf (float y, float *sinp, float *cosp)
+{
+ double x = y;
+ double s;
+ int n;
+ sincos_t *p = &sincosf_table[0];
+
+ if (abstop12 (y) < abstop12 (pio4))
+ {
+ double x2 = x * x;
+
+ if (unlikely (abstop12 (y) < abstop12 (0x1p-12f)))
+ {
+ if (unlikely (abstop12 (y) < abstop12 (0x1p-126f)))
+ /* Force underflow for tiny y. */
+ force_eval_float (x2);
+ *sinp = y;
+ *cosp = 1.0f;
+ return;
+ }
+
+ sincosf_poly (x, x2, p, 0, sinp, cosp);
+ }
+ else if (abstop12 (y) < abstop12 (120.0f))
+ {
+ x = reduce_fast (x, p, &n);
+
+ /* Setup the signs for sin and cos. */
+ s = p->sign[n & 3];
+
+ if (n & 2)
+ p = &sincosf_table[1];
+
+ sincosf_poly (x * s, x * x, p, n, sinp, cosp);
+ }
+ else if (likely (abstop12 (y) < abstop12 (INFINITY)))
+ {
+ uint32_t xi = asuint (y);
+ int sign = xi >> 31;
+
+ x = reduce_large (xi, &n);
+
+ /* Setup signs for sin and cos - include original sign. */
+ s = p->sign[(n + sign) & 3];
+
+ if ((n + sign) & 2)
+ p = &sincosf_table[1];
+
+ sincosf_poly (x * s, x * x, p, n, sinp, cosp);
+ }
+ else
+ {
+ /* Return NaN if Inf or NaN for both sin and cos. */
+ *sinp = *cosp = y - y;
+#if WANT_ERRNO
+ /* Needed to set errno for +-Inf, the add is a hack to work
+ around a gcc register allocation issue: just passing y
+ affects code generation in the fast path. */
+ __math_invalidf (y + y);
+#endif
+ }
+}
+
+#endif
diff --git a/newlib/libm/common/sf_sincos_data.c b/newlib/libm/common/sf_sincos_data.c
new file mode 100644
index 0000000000000000000000000000000000000000..cecc1f393268b74bc9101142a2915688f352c87c
--- /dev/null
+++ b/newlib/libm/common/sf_sincos_data.c
@@ -0,0 +1,87 @@
+/* Data definitions for sinf, cosf and sincosf.
+ Copyright (c) 2018 Arm Ltd. All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+ 1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+ 2. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+ 3. The name of the company may not be used to endorse or promote
+ products derived from this software without specific prior written
+ permission.
+
+ THIS SOFTWARE IS PROVIDED BY ARM LTD ``AS IS AND ANY EXPRESS OR IMPLIED
+ WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+ MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ IN NO EVENT SHALL ARM LTD BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+ SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
+ TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */
+
+#include "fdlibm.h"
+#if !__OBSOLETE_MATH
+
+#include <stdint.h>
+#include <math.h>
+#include "math_config.h"
+#include "sf_sincos.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 TOINT_INTRINSICS
+ 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 TOINT_INTRINSICS
+ 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
diff --git a/newlib/libm/math/sf_cos.c b/newlib/libm/math/sf_cos.c
index 4c0a9a53569ada72f1449517297d116af32534aa..8cb0eb0087a36e5763db869aa91ee6731a2ed7d9 100644
--- a/newlib/libm/math/sf_cos.c
+++ b/newlib/libm/math/sf_cos.c
@@ -14,6 +14,7 @@
*/
#include "fdlibm.h"
+#if __OBSOLETE_MATH
#ifdef __STDC__
static const float one=1.0;
@@ -66,3 +67,4 @@ static float one=1.0;
}
#endif /* defined(_DOUBLE_IS_32BITS) */
+#endif /* __OBSOLETE_MATH */
diff --git a/newlib/libm/math/sf_sin.c b/newlib/libm/math/sf_sin.c
index da81845d98ef125549baf51e9e1cb00784d98499..c3ec6e2e40bd763828d103ba6fc516711532729e 100644
--- a/newlib/libm/math/sf_sin.c
+++ b/newlib/libm/math/sf_sin.c
@@ -14,6 +14,7 @@
*/
#include "fdlibm.h"
+#if __OBSOLETE_MATH
#ifdef __STDC__
float sinf(float x)
@@ -60,3 +61,4 @@
}
#endif /* defined(_DOUBLE_IS_32BITS) */
+#endif /* __OBSOLETE_MATH */
diff --git a/newlib/libm/math/wf_sincos.c b/newlib/libm/math/wf_sincos.c
index 477c604019f750b437aa6b0d3e8c8733e7634a0d..69eb922c982707b94f3eaa2da509d35dbb16faae 100644
--- a/newlib/libm/math/wf_sincos.c
+++ b/newlib/libm/math/wf_sincos.c
@@ -1,6 +1,8 @@
/* sincos -- currently no more efficient than two separate calls to
sin and cos. */
#include "fdlibm.h"
+#if __OBSOLETE_MATH
+
#include <errno.h>
#ifdef __STDC__
@@ -31,3 +33,4 @@
*cosx = cosf((float) x);
}
#endif /* defined(_DOUBLE_IS_32BITS) */
+#endif /* __OBSOLETE_MATH */