Inverse sine function.
Method
Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
we approximate asin(x) on [0,0.5] by
asin(x) = x + x*x^2*R(x^2)
where
R(x^2) is a rational approximation of (asin(x)-x)/x^3
and its remez error is bounded by
|(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)For x in [0.5,1]
asin(x) = pi/2-2*asin(sqrt((1-x)/2))
Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
then for x>0.98
asin(x) = pi/2 - 2*(s+s*z*R(z))
= pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)For x<=0.98, let pio4_hi = pio2_hi/2, then
f = hi part of s;
c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z)
and
asin(x) = pi/2 - 2*(s+s*z*R(z))
= pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
= pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))Special cases
if x is NaN, return NaN
if |x|>1, return NaN
Definition at line 100 of file asin.c.
References fabs(), huge, one, pio2_hi, pio2_lo, pio4_hi, pS0, pS1, pS2, pS3, pS4, pS5, qS1, qS2, qS3, qS4, and sqrt().
Referenced by acrd(), acsc(), acvc(), aexc(), and RealNumber::ArcSine().
102 double t, w, p, q, c, r, s;
105 ix = hx & 0x7fffffff;
106 if (ix >= 0x3ff00000)
110 if (((ix - 0x3ff00000) | lx) == 0)
115 else if (ix < 0x3fe00000)
144 if (ix >= 0x3FEF3333)
153 c = (t - w * w) / (s + w);
155 p = 2.0 * s * r - (
pio2_lo - 2.0 * c);
#define GET_HIGH_WORD(i, d)
Get the more significant 32 bit int from a double.
static const double pio2_hi
#define GET_LOW_WORD(i, d)
Get the less significant 32 bit int from a double.
static const double pio4_hi
static const double pio2_lo
#define SET_LOW_WORD(d, v)
Set the less significant 32 bits of a double from an int.
double sqrt(double x)
Square root function.
double fabs(double x)
Returns the absolute value of x.
