Back to index

glibc  2.9
e_asinl.c
Go to the documentation of this file.
00001 /*
00002  * ====================================================
00003  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00004  *
00005  * Developed at SunPro, a Sun Microsystems, Inc. business.
00006  * Permission to use, copy, modify, and distribute this
00007  * software is freely granted, provided that this notice
00008  * is preserved.
00009  * ====================================================
00010  */
00011 
00012 /*
00013   Long double expansions are
00014   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
00015   and are incorporated herein by permission of the author.  The author 
00016   reserves the right to distribute this material elsewhere under different
00017   copying permissions.  These modifications are distributed here under 
00018   the following terms:
00019 
00020     This library is free software; you can redistribute it and/or
00021     modify it under the terms of the GNU Lesser General Public
00022     License as published by the Free Software Foundation; either
00023     version 2.1 of the License, or (at your option) any later version.
00024 
00025     This library is distributed in the hope that it will be useful,
00026     but WITHOUT ANY WARRANTY; without even the implied warranty of
00027     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00028     Lesser General Public License for more details.
00029 
00030     You should have received a copy of the GNU Lesser General Public
00031     License along with this library; if not, write to the Free Software
00032     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00033 
00034 /* __ieee754_asin(x)
00035  * Method :
00036  *     Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
00037  *     we approximate asin(x) on [0,0.5] by
00038  *            asin(x) = x + x*x^2*R(x^2)
00039  *
00040  *     For x in [0.5,1]
00041  *            asin(x) = pi/2-2*asin(sqrt((1-x)/2))
00042  *     Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
00043  *     then for x>0.98
00044  *            asin(x) = pi/2 - 2*(s+s*z*R(z))
00045  *                   = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
00046  *     For x<=0.98, let pio4_hi = pio2_hi/2, then
00047  *            f = hi part of s;
00048  *            c = sqrt(z) - f = (z-f*f)/(s+f)    ...f+c=sqrt(z)
00049  *     and
00050  *            asin(x) = pi/2 - 2*(s+s*z*R(z))
00051  *                   = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
00052  *                   = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
00053  *
00054  * Special cases:
00055  *     if x is NaN, return x itself;
00056  *     if |x|>1, return NaN with invalid signal.
00057  *
00058  */
00059 
00060 
00061 #include "math.h"
00062 #include "math_private.h"
00063 
00064 #ifdef __STDC__
00065 static const long double
00066 #else
00067 static long double
00068 #endif
00069   one = 1.0L,
00070   huge = 1.0e+4932L,
00071  pio2_hi = 1.5707963267948966192021943710788178805159986950457096099853515625L,
00072   pio2_lo = 2.9127320560933561582586004641843300502121E-20L,
00073   pio4_hi = 7.8539816339744830960109718553940894025800E-1L,
00074 
00075        /* coefficient for R(x^2) */
00076 
00077   /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
00078      0 <= x <= 0.5
00079      peak relative error 1.9e-21  */
00080   pS0 =  -1.008714657938491626019651170502036851607E1L,
00081   pS1 =   2.331460313214179572063441834101394865259E1L,
00082   pS2 =  -1.863169762159016144159202387315381830227E1L,
00083   pS3 =   5.930399351579141771077475766877674661747E0L,
00084   pS4 =  -6.121291917696920296944056882932695185001E-1L,
00085   pS5 =   3.776934006243367487161248678019350338383E-3L,
00086 
00087   qS0 =  -6.052287947630949712886794360635592886517E1L,
00088   qS1 =   1.671229145571899593737596543114258558503E2L,
00089   qS2 =  -1.707840117062586426144397688315411324388E2L,
00090   qS3 =   7.870295154902110425886636075950077640623E1L,
00091   qS4 =  -1.568433562487314651121702982333303458814E1L;
00092     /* 1.000000000000000000000000000000000000000E0 */
00093 
00094 #ifdef __STDC__
00095 long double
00096 __ieee754_asinl (long double x)
00097 #else
00098 double
00099 __ieee754_asinl (x)
00100      long double x;
00101 #endif
00102 {
00103   long double t, w, p, q, c, r, s;
00104   int32_t ix;
00105   u_int32_t se, i0, i1, k;
00106 
00107   GET_LDOUBLE_WORDS (se, i0, i1, x);
00108   ix = se & 0x7fff;
00109   ix = (ix << 16) | (i0 >> 16);
00110   if (ix >= 0x3fff8000)
00111     {                       /* |x|>= 1 */
00112       if (ix == 0x3fff8000 && ((i0 - 0x80000000) | i1) == 0)
00113        /* asin(1)=+-pi/2 with inexact */
00114        return x * pio2_hi + x * pio2_lo;
00115       return (x - x) / (x - x);    /* asin(|x|>1) is NaN */
00116     }
00117   else if (ix < 0x3ffe8000)
00118     {                       /* |x|<0.5 */
00119       if (ix < 0x3fde8000)
00120        {                    /* if |x| < 2**-33 */
00121          if (huge + x > one)
00122            return x;        /* return x with inexact if x!=0 */
00123        }
00124       else
00125        {
00126          t = x * x;
00127          p =
00128            t * (pS0 +
00129                t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
00130          q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t))));
00131          w = p / q;
00132          return x + x * w;
00133        }
00134     }
00135   /* 1> |x|>= 0.5 */
00136   w = one - fabsl (x);
00137   t = w * 0.5;
00138   p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
00139   q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t))));
00140   s = __ieee754_sqrtl (t);
00141   if (ix >= 0x3ffef999)
00142     {                       /* if |x| > 0.975 */
00143       w = p / q;
00144       t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
00145     }
00146   else
00147     {
00148       GET_LDOUBLE_WORDS (k, i0, i1, s);
00149       i1 = 0;
00150       SET_LDOUBLE_WORDS (w,k,i0,i1);
00151       c = (t - w * w) / (s + w);
00152       r = p / q;
00153       p = 2.0 * s * r - (pio2_lo - 2.0 * c);
00154       q = pio4_hi - 2.0 * w;
00155       t = pio4_hi - (p - q);
00156     }
00157   if ((se & 0x8000) == 0)
00158     return t;
00159   else
00160     return -t;
00161 }