PlaneteCasio
/
libc
1
0
Fork 0
Code Issues 6 Pull Requests Releases Wiki Activity

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.

--
This commit is contained in:
Wilco Dijkstra 2018-06-20 12:07:22 +00:00 committed by Corinna Vinschen
parent cfe8c6c504
commit 3baadb9912
11 changed files with 667 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
newlib/libm/common/Makefile.am
View File

@ -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

34
newlib/libm/common/Makefile.in
View File

@ -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

88
newlib/libm/common/cosf.c Normal file
View File

@ -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

84
newlib/libm/common/math_config.h
View File

@ -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);

104
newlib/libm/common/sincosf.c Normal file
View File

@ -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

172
newlib/libm/common/sincosf.h Normal file
View File

@ -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;
}

87
newlib/libm/common/sincosf_data.c Normal file
View File

@ -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

92
newlib/libm/common/sinf.c Normal file
View File

@ -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

2
newlib/libm/math/sf_cos.c
View File

@ -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 */

2
newlib/libm/math/sf_sin.c
View File

@ -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 */

3
newlib/libm/math/wf_sincos.c
View File

@ -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 */