Revert "Improve performance of sinf/cosf/sincosf"

This reverts commit fca80a9d1b.

Accidentally pushed a preliminary version
This commit is contained in:
Corinna Vinschen 2018-06-21 09:36:11 +02:00
parent b7d9d27b0e
commit cfe8c6c504
11 changed files with 6 additions and 667 deletions

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 \
sinf.c cosf.c sincosf.c sincosf_data.c math_errf.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,7 +38,6 @@ 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
@ -53,7 +52,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

View File

@ -55,7 +55,6 @@ 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
@ -115,8 +114,6 @@ 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) \
@ -172,8 +169,7 @@ 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 sinf.lo cosf.lo sincosf.lo sincosf_data.lo \
math_errf.lo
sf_pow.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 \
@ -365,7 +361,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 \
sinf.c cosf.c sincosf.c sincosf_data.c math_errf.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 \
@ -378,7 +374,6 @@ 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)
@ -386,6 +381,7 @@ lib_a_CFLAGS = -fbuiltin -fno-math-errno $(am__append_3)
@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
@ -948,30 +944,6 @@ 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

View File

@ -1,88 +0,0 @@
/* 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

View File

@ -1,5 +1,5 @@
/* Configuration for math routines.
Copyright (c) 2017-2018 Arm Ltd. All rights reserved.
Copyright (c) 2017 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,44 +43,6 @@
# 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
@ -147,56 +109,12 @@ 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);

View File

@ -1,104 +0,0 @@
/* 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

View File

@ -1,172 +0,0 @@
/* 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;
}

View File

@ -1,87 +0,0 @@
/* 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

View File

@ -1,92 +0,0 @@
/* 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

View File

@ -14,7 +14,6 @@
*/
#include "fdlibm.h"
#if __OBSOLETE_MATH
#ifdef __STDC__
static const float one=1.0;
@ -67,4 +66,3 @@ static float one=1.0;
}
#endif /* defined(_DOUBLE_IS_32BITS) */
#endif /* __OBSOLETE_MATH */

View File

@ -14,7 +14,6 @@
*/
#include "fdlibm.h"
#if __OBSOLETE_MATH
#ifdef __STDC__
float sinf(float x)
@ -61,4 +60,3 @@
}
#endif /* defined(_DOUBLE_IS_32BITS) */
#endif /* __OBSOLETE_MATH */

View File

@ -1,8 +1,6 @@
/* sincos -- currently no more efficient than two separate calls to
sin and cos. */
#include "fdlibm.h"
#if __OBSOLETE_MATH
#include <errno.h>
#ifdef __STDC__
@ -33,4 +31,3 @@
*cosx = cosf((float) x);
}
#endif /* defined(_DOUBLE_IS_32BITS) */
#endif /* __OBSOLETE_MATH */