public inbox for newlib-cvs@sourceware.org
help / color / mirror / Atom feed
* [newlib-cygwin] Improve performance of sinf/cosf/sincosf
@ 2018-06-21 7:35 Corinna Vinschen
0 siblings, 0 replies; 2+ messages in thread
From: Corinna Vinschen @ 2018-06-21 7:35 UTC (permalink / raw)
To: newlib-cvs
https://sourceware.org/git/gitweb.cgi?p=newlib-cygwin.git;h=fca80a9d1b3fa6620cdaccec6b726eef1a6530a1
commit fca80a9d1b3fa6620cdaccec6b726eef1a6530a1
Author: Wilco Dijkstra <Wilco.Dijkstra@arm.com>
Date: Mon Jun 18 17:59:37 2018 +0000
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.
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-06-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:
---
newlib/libm/common/Makefile.am | 5 +-
newlib/libm/common/Makefile.in | 34 +++++++-
newlib/libm/common/cosf.c | 88 +++++++++++++++++++
newlib/libm/common/math_config.h | 84 ++++++++++++++++++-
newlib/libm/common/sincosf.c | 104 +++++++++++++++++++++++
newlib/libm/common/sincosf.h | 172 ++++++++++++++++++++++++++++++++++++++
newlib/libm/common/sincosf_data.c | 87 +++++++++++++++++++
newlib/libm/common/sinf.c | 92 ++++++++++++++++++++
newlib/libm/math/sf_cos.c | 2 +
newlib/libm/math/sf_sin.c | 2 +
newlib/libm/math/wf_sincos.c | 3 +
11 files changed, 667 insertions(+), 6 deletions(-)
diff --git a/newlib/libm/common/Makefile.am b/newlib/libm/common/Makefile.am
index 3e2620c..97b2431 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
+ sinf.c cosf.c sincosf.c sincosf_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 6a0ba73..c944ac7 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-sinf.$(OBJEXT) lib_a-cosf.$(OBJEXT) \
+ lib_a-sincosf.$(OBJEXT) lib_a-sincosf_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 sinf.lo cosf.lo sincosf.lo sincosf_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
+ sinf.c cosf.c sincosf.c sincosf_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-sinf.o: sinf.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sinf.o `test -f 'sinf.c' || echo '$(srcdir)/'`sinf.c
+
+lib_a-sinf.obj: sinf.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sinf.obj `if test -f 'sinf.c'; then $(CYGPATH_W) 'sinf.c'; else $(CYGPATH_W) '$(srcdir)/sinf.c'; fi`
+
+lib_a-cosf.o: cosf.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-cosf.o `test -f 'cosf.c' || echo '$(srcdir)/'`cosf.c
+
+lib_a-cosf.obj: cosf.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-cosf.obj `if test -f 'cosf.c'; then $(CYGPATH_W) 'cosf.c'; else $(CYGPATH_W) '$(srcdir)/cosf.c'; fi`
+
+lib_a-sincosf.o: sincosf.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sincosf.o `test -f 'sincosf.c' || echo '$(srcdir)/'`sincosf.c
+
+lib_a-sincosf.obj: sincosf.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sincosf.obj `if test -f 'sincosf.c'; then $(CYGPATH_W) 'sincosf.c'; else $(CYGPATH_W) '$(srcdir)/sincosf.c'; fi`
+
+lib_a-sincosf_data.o: sincosf_data.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sincosf_data.o `test -f 'sincosf_data.c' || echo '$(srcdir)/'`sincosf_data.c
+
+lib_a-sincosf_data.obj: sincosf_data.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sincosf_data.obj `if test -f 'sincosf_data.c'; then $(CYGPATH_W) 'sincosf_data.c'; else $(CYGPATH_W) '$(srcdir)/sincosf_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/cosf.c b/newlib/libm/common/cosf.c
new file mode 100644
index 0000000..aac0a9a
--- /dev/null
+++ b/newlib/libm/common/cosf.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 "sincosf.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/math_config.h b/newlib/libm/common/math_config.h
index c566784..6d87376 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/sincosf.c b/newlib/libm/common/sincosf.c
new file mode 100644
index 0000000..85f1264
--- /dev/null
+++ b/newlib/libm/common/sincosf.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 "sincosf.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/sincosf.h b/newlib/libm/common/sincosf.h
new file mode 100644
index 0000000..2f3b2f8
--- /dev/null
+++ b/newlib/libm/common/sincosf.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 = ((int)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/sincosf_data.c b/newlib/libm/common/sincosf_data.c
new file mode 100644
index 0000000..069ee7b
--- /dev/null
+++ b/newlib/libm/common/sincosf_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 "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 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/common/sinf.c b/newlib/libm/common/sinf.c
new file mode 100644
index 0000000..c1baf42
--- /dev/null
+++ b/newlib/libm/common/sinf.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 "sincosf.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/math/sf_cos.c b/newlib/libm/math/sf_cos.c
index 4c0a9a5..8cb0eb0 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 da81845..c3ec6e2 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 477c604..69eb922 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 */
^ permalink raw reply [flat|nested] 2+ messages in thread
* [newlib-cygwin] Improve performance of sinf/cosf/sincosf
@ 2018-06-21 7:39 Corinna Vinschen
0 siblings, 0 replies; 2+ messages in thread
From: Corinna Vinschen @ 2018-06-21 7:39 UTC (permalink / raw)
To: newlib-cvs
https://sourceware.org/git/gitweb.cgi?p=newlib-cygwin.git;h=3baadb9912c11bbd0e8067b38ae2e2a9e206cbb7
commit 3baadb9912c11bbd0e8067b38ae2e2a9e206cbb7
Author: Wilco Dijkstra <Wilco.Dijkstra@arm.com>
Date: Wed Jun 20 12:07:22 2018 +0000
Improve performance of sinf/cosf/sincosf
Here is the correct patch with both filenames and int cast fixed:
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:
---
newlib/libm/common/Makefile.am | 5 +-
newlib/libm/common/Makefile.in | 34 +++++++-
newlib/libm/common/cosf.c | 88 +++++++++++++++++++
newlib/libm/common/math_config.h | 84 ++++++++++++++++++-
newlib/libm/common/sincosf.c | 104 +++++++++++++++++++++++
newlib/libm/common/sincosf.h | 172 ++++++++++++++++++++++++++++++++++++++
newlib/libm/common/sincosf_data.c | 87 +++++++++++++++++++
newlib/libm/common/sinf.c | 92 ++++++++++++++++++++
newlib/libm/math/sf_cos.c | 2 +
newlib/libm/math/sf_sin.c | 2 +
newlib/libm/math/wf_sincos.c | 3 +
11 files changed, 667 insertions(+), 6 deletions(-)
diff --git a/newlib/libm/common/Makefile.am b/newlib/libm/common/Makefile.am
index 3e2620c..97b2431 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
+ sinf.c cosf.c sincosf.c sincosf_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 6a0ba73..c944ac7 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-sinf.$(OBJEXT) lib_a-cosf.$(OBJEXT) \
+ lib_a-sincosf.$(OBJEXT) lib_a-sincosf_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 sinf.lo cosf.lo sincosf.lo sincosf_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
+ sinf.c cosf.c sincosf.c sincosf_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-sinf.o: sinf.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sinf.o `test -f 'sinf.c' || echo '$(srcdir)/'`sinf.c
+
+lib_a-sinf.obj: sinf.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sinf.obj `if test -f 'sinf.c'; then $(CYGPATH_W) 'sinf.c'; else $(CYGPATH_W) '$(srcdir)/sinf.c'; fi`
+
+lib_a-cosf.o: cosf.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-cosf.o `test -f 'cosf.c' || echo '$(srcdir)/'`cosf.c
+
+lib_a-cosf.obj: cosf.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-cosf.obj `if test -f 'cosf.c'; then $(CYGPATH_W) 'cosf.c'; else $(CYGPATH_W) '$(srcdir)/cosf.c'; fi`
+
+lib_a-sincosf.o: sincosf.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sincosf.o `test -f 'sincosf.c' || echo '$(srcdir)/'`sincosf.c
+
+lib_a-sincosf.obj: sincosf.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sincosf.obj `if test -f 'sincosf.c'; then $(CYGPATH_W) 'sincosf.c'; else $(CYGPATH_W) '$(srcdir)/sincosf.c'; fi`
+
+lib_a-sincosf_data.o: sincosf_data.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sincosf_data.o `test -f 'sincosf_data.c' || echo '$(srcdir)/'`sincosf_data.c
+
+lib_a-sincosf_data.obj: sincosf_data.c
+ $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(lib_a_CFLAGS) $(CFLAGS) -c -o lib_a-sincosf_data.obj `if test -f 'sincosf_data.c'; then $(CYGPATH_W) 'sincosf_data.c'; else $(CYGPATH_W) '$(srcdir)/sincosf_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/cosf.c b/newlib/libm/common/cosf.c
new file mode 100644
index 0000000..aac0a9a
--- /dev/null
+++ b/newlib/libm/common/cosf.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 "sincosf.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/math_config.h b/newlib/libm/common/math_config.h
index c566784..6d87376 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/sincosf.c b/newlib/libm/common/sincosf.c
new file mode 100644
index 0000000..85f1264
--- /dev/null
+++ b/newlib/libm/common/sincosf.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 "sincosf.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/sincosf.h b/newlib/libm/common/sincosf.h
new file mode 100644
index 0000000..955a733
--- /dev/null
+++ b/newlib/libm/common/sincosf.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/sincosf_data.c b/newlib/libm/common/sincosf_data.c
new file mode 100644
index 0000000..069ee7b
--- /dev/null
+++ b/newlib/libm/common/sincosf_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 "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 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/common/sinf.c b/newlib/libm/common/sinf.c
new file mode 100644
index 0000000..c1baf42
--- /dev/null
+++ b/newlib/libm/common/sinf.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 "sincosf.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/math/sf_cos.c b/newlib/libm/math/sf_cos.c
index 4c0a9a5..8cb0eb0 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 da81845..c3ec6e2 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 477c604..69eb922 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 */
^ permalink raw reply [flat|nested] 2+ messages in thread
end of thread, other threads:[~2018-06-21 7:39 UTC | newest]
Thread overview: 2+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2018-06-21 7:35 [newlib-cygwin] Improve performance of sinf/cosf/sincosf Corinna Vinschen
2018-06-21 7:39 Corinna Vinschen
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).