sysdeps/ieee754/ldbl-96/e_asinl.c

   1 /*
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10  */
  11
  12 /*
  13   Long double expansions are
  14   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
  15   and are incorporated herein by permission of the author.  The author
  16   reserves the right to distribute this material elsewhere under different
  17   copying permissions.  These modifications are distributed here under
  18   the following terms:
  19
  20     This library is free software; you can redistribute it and/or
  21     modify it under the terms of the GNU Lesser General Public
  22     License as published by the Free Software Foundation; either
  23     version 2.1 of the License, or (at your option) any later version.
  24
  25     This library is distributed in the hope that it will be useful,
  26     but WITHOUT ANY WARRANTY; without even the implied warranty of
  27     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  28     Lesser General Public License for more details.
  29
  30     You should have received a copy of the GNU Lesser General Public
  31     License along with this library; if not, write to the Free Software
  32     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
  33
  34 /* __ieee754_asin(x)
  35  * Method :
  36  *      Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
  37  *      we approximate asin(x) on [0,0.5] by
  38  *              asin(x) = x + x*x^2*R(x^2)
  39  *
  40  *      For x in [0.5,1]
  41  *              asin(x) = pi/2-2*asin(sqrt((1-x)/2))
  42  *      Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
  43  *      then for x>0.98
  44  *              asin(x) = pi/2 - 2*(s+s*z*R(z))
  45  *                      = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
  46  *      For x<=0.98, let pio4_hi = pio2_hi/2, then
  47  *              f = hi part of s;
  48  *              c = sqrt(z) - f = (z-f*f)/(s+f)         ...f+c=sqrt(z)
  49  *      and
  50  *              asin(x) = pi/2 - 2*(s+s*z*R(z))
  51  *                      = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
  52  *                      = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
  53  *
  54  * Special cases:
  55  *      if x is NaN, return x itself;
  56  *      if |x|>1, return NaN with invalid signal.
  57  *
  58  */
  59
  60
  61 #include "math.h"
  62 #include "math_private.h"
  63
  64 #ifdef __STDC__
  65 static const long double
  66 #else
  67 static long double
  68 #endif
  69   one = 1.0L,
  70   huge = 1.0e+4932L,
  71  pio2_hi = 1.5707963267948966192021943710788178805159986950457096099853515625L,
  72   pio2_lo = 2.9127320560933561582586004641843300502121E-20L,
  73   pio4_hi = 7.8539816339744830960109718553940894025800E-1L,
  74
  75         /* coefficient for R(x^2) */
  76
  77   /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
  78      0 <= x <= 0.5
  79      peak relative error 1.9e-21  */
  80   pS0 =  -1.008714657938491626019651170502036851607E1L,
  81   pS1 =   2.331460313214179572063441834101394865259E1L,
  82   pS2 =  -1.863169762159016144159202387315381830227E1L,
  83   pS3 =   5.930399351579141771077475766877674661747E0L,
  84   pS4 =  -6.121291917696920296944056882932695185001E-1L,
  85   pS5 =   3.776934006243367487161248678019350338383E-3L,
  86
  87   qS0 =  -6.052287947630949712886794360635592886517E1L,
  88   qS1 =   1.671229145571899593737596543114258558503E2L,
  89   qS2 =  -1.707840117062586426144397688315411324388E2L,
  90   qS3 =   7.870295154902110425886636075950077640623E1L,
  91   qS4 =  -1.568433562487314651121702982333303458814E1L;
  92     /* 1.000000000000000000000000000000000000000E0 */
  93
  94 #ifdef __STDC__
  95 long double
  96 __ieee754_asinl (long double x)
  97 #else
  98 double
  99 __ieee754_asinl (x)
 100      long double x;
 101 #endif
 102 {
 103   long double t, w, p, q, c, r, s;
 104   int32_t ix;
 105   u_int32_t se, i0, i1, k;
 106
 107   GET_LDOUBLE_WORDS (se, i0, i1, x);
 108   ix = se & 0x7fff;
 109   ix = (ix << 16) | (i0 >> 16);
 110   if (ix >= 0x3fff8000)
 111     {                           /* |x|>= 1 */
 112       if (ix == 0x3fff8000 && ((i0 - 0x80000000) | i1) == 0)
 113         /* asin(1)=+-pi/2 with inexact */
 114         return x * pio2_hi + x * pio2_lo;
 115       return (x - x) / (x - x); /* asin(|x|>1) is NaN */
 116     }
 117   else if (ix < 0x3ffe8000)
 118     {                           /* |x|<0.5 */
 119       if (ix < 0x3fde8000)
 120         {                       /* if |x| < 2**-33 */
 121           if (huge + x > one)
 122             return x;           /* return x with inexact if x!=0 */
 123         }
 124       else
 125         {
 126           t = x * x;
 127           p =
 128             t * (pS0 +
 129                  t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
 130           q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t))));
 131           w = p / q;
 132           return x + x * w;
 133         }
 134     }
 135   /* 1> |x|>= 0.5 */
 136   w = one - fabsl (x);
 137   t = w * 0.5;
 138   p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
 139   q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t))));
 140   s = __ieee754_sqrtl (t);
 141   if (ix >= 0x3ffef999)
 142     {                           /* if |x| > 0.975 */
 143       w = p / q;
 144       t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
 145     }
 146   else
 147     {
 148       GET_LDOUBLE_WORDS (k, i0, i1, s);
 149       i1 = 0;
 150       SET_LDOUBLE_WORDS (w,k,i0,i1);
 151       c = (t - w * w) / (s + w);
 152       r = p / q;
 153       p = 2.0 * s * r - (pio2_lo - 2.0 * c);
 154       q = pio4_hi - 2.0 * w;
 155       t = pio4_hi - (p - q);
 156     }
 157   if ((se & 0x8000) == 0)
 158     return t;
 159   else
 160     return -t;
 161 }