/*
|
* rredf.h - trigonometric range reduction function written new for RVCT 4.1
|
*
|
* Copyright (c) 2009-2018, Arm Limited.
|
* SPDX-License-Identifier: MIT
|
*/
|
|
/*
|
* This header file defines an inline function which all three of
|
* the single-precision trig functions (sinf, cosf, tanf) should use
|
* to perform range reduction. The inline function handles the
|
* quickest and most common cases inline, before handing off to an
|
* out-of-line function defined in rredf.c for everything else. Thus
|
* a reasonable compromise is struck between speed and space. (I
|
* hope.) In particular, this approach avoids a function call
|
* overhead in the common case.
|
*/
|
|
#ifndef _included_rredf_h
|
#define _included_rredf_h
|
|
#include "math_private.h"
|
|
#ifdef __cplusplus
|
extern "C" {
|
#endif /* __cplusplus */
|
|
extern float __mathlib_rredf2(float x, int *q, unsigned k);
|
|
/*
|
* Semantics of the function:
|
* - x is the single-precision input value provided by the user
|
* - the return value is in the range [-pi/4,pi/4], and is equal
|
* (within reasonable accuracy bounds) to x minus n*pi/2 for some
|
* integer n. (FIXME: perhaps some slippage on the output
|
* interval is acceptable, requiring more range from the
|
* following polynomial approximations but permitting more
|
* approximate rred decisions?)
|
* - *q is set to a positive value whose low two bits match those
|
* of n. Alternatively, it comes back as -1 indicating that the
|
* input value was trivial in some way (infinity, NaN, or so
|
* small that we can safely return sin(x)=tan(x)=x,cos(x)=1).
|
*/
|
static __inline float __mathlib_rredf(float x, int *q)
|
{
|
/*
|
* First, extract the bit pattern of x as an integer, so that we
|
* can repeatedly compare things to it without multiple
|
* overheads in retrieving comparison results from the VFP.
|
*/
|
unsigned k = fai(x);
|
|
/*
|
* Deal immediately with the simplest possible case, in which x
|
* is already within the interval [-pi/4,pi/4]. This also
|
* identifies the subcase of ludicrously small x.
|
*/
|
if ((k << 1) < (0x3f490fdb << 1)) {
|
if ((k << 1) < (0x39800000 << 1))
|
*q = -1;
|
else
|
*q = 0;
|
return x;
|
}
|
|
/*
|
* The next plan is to multiply x by 2/pi and convert to an
|
* integer, which gives us n; then we subtract n*pi/2 from x to
|
* get our output value.
|
*
|
* By representing pi/2 in that final step by a prec-and-a-half
|
* approximation, we can arrange good accuracy for n strictly
|
* less than 2^13 (so that an FP representation of n has twelve
|
* zero bits at the bottom). So our threshold for this strategy
|
* is 2^13 * pi/2 - pi/4, otherwise known as 8191.75 * pi/2 or
|
* 4095.875*pi. (Or, for those perverse people interested in
|
* actual numbers rather than multiples of pi/2, about 12867.5.)
|
*/
|
if (__builtin_expect((k & 0x7fffffff) < 0x46490e49, 1)) {
|
float nf = 0.636619772367581343f * x;
|
/*
|
* The difference between that single-precision constant and
|
* the real 2/pi is about 2.568e-8. Hence the product nf has a
|
* potential error of 2.568e-8|x| even before rounding; since
|
* |x| < 4096 pi, that gives us an error bound of about
|
* 0.0003305.
|
*
|
* nf is then rounded to single precision, with a max error of
|
* 1/2 ULP, and since nf goes up to just under 8192, half a
|
* ULP could be as big as 2^-12 ~= 0.0002441.
|
*
|
* So by the time we convert nf to an integer, it could be off
|
* by that much, causing the wrong integer to be selected, and
|
* causing us to return a value a little bit outside the
|
* theoretical [-pi/4,+pi/4] output interval.
|
*
|
* How much outside? Well, we subtract nf*pi/2 from x, so the
|
* error bounds above have be be multiplied by pi/2. And if
|
* both of the above sources of error suffer their worst cases
|
* at once, then the very largest value we could return is
|
* obtained by adding that lot to the interval bound pi/4 to
|
* get
|
*
|
* pi/4 + ((2/pi - 0f_3f22f983)*4096*pi + 2^-12) * pi/2
|
*
|
* which comes to 0f_3f494b02. (Compare 0f_3f490fdb = pi/4.)
|
*
|
* So callers of this range reducer should be prepared to
|
* handle numbers up to that large.
|
*/
|
#ifdef __TARGET_FPU_SOFTVFP
|
nf = _frnd(nf);
|
#else
|
if (k & 0x80000000)
|
nf = (nf - 8388608.0f) + 8388608.0f;
|
else
|
nf = (nf + 8388608.0f) - 8388608.0f; /* round to _nearest_ integer. FIXME: use some sort of frnd in softfp */
|
#endif
|
*q = 3 & (int)nf;
|
#if 0
|
/*
|
* FIXME: now I need a bunch of special cases to avoid
|
* having to do the full four-word reduction every time.
|
* Also, adjust the comment at the top of this section!
|
*/
|
if (__builtin_expect((k & 0x7fffffff) < 0x46490e49, 1))
|
return ((x - nf * 0x1.92p+0F) - nf * 0x1.fb4p-12F) - nf * 0x1.4442d2p-24F;
|
else
|
#endif
|
return ((x - nf * 0x1.92p+0F) - nf * 0x1.fb4p-12F) - nf * 0x1.444p-24F - nf * 0x1.68c234p-39F;
|
}
|
|
/*
|
* That's enough to do in-line; if we're still playing, hand off
|
* to the out-of-line main range reducer.
|
*/
|
return __mathlib_rredf2(x, q, k);
|
}
|
|
#ifdef __cplusplus
|
} /* end of extern "C" */
|
#endif /* __cplusplus */
|
|
#endif /* included */
|
|
/* end of rredf.h */
|
