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, see
  32     <http://www.gnu.org/licenses/>.  */
  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 <float.h>
  62 #include <math.h>
  63 #include <math_private.h>
  64
  65 static const long double
  66   one = 1.0L,
  67   huge = 1.0e+4932L,
  68   pio2_hi = 0x1.921fb54442d1846ap+0L, /* pi/2 rounded to nearest to 64
  69                                          bits.  */
  70   pio2_lo = -0x7.6733ae8fe47c65d8p-68L, /* pi/2 - pio2_hi rounded to
  71                                            nearest to 64 bits.  */
  72   pio4_hi = 0xc.90fdaa22168c235p-4L, /* pi/4 rounded to nearest to 64
  73                                         bits.  */
  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 long double
  95 __ieee754_asinl (long double x)
  96 {
  97   long double t, w, p, q, c, r, s;
  98   int32_t ix;
  99   u_int32_t se, i0, i1, k;
 100
 101   GET_LDOUBLE_WORDS (se, i0, i1, x);
 102   ix = se & 0x7fff;
 103   ix = (ix << 16) | (i0 >> 16);
 104   if (ix >= 0x3fff8000)
 105     {                           /* |x|>= 1 */
 106       if (ix == 0x3fff8000 && ((i0 - 0x80000000) | i1) == 0)
 107         /* asin(1)=+-pi/2 with inexact */
 108         return x * pio2_hi + x * pio2_lo;
 109       return (x - x) / (x - x); /* asin(|x|>1) is NaN */
 110     }
 111   else if (ix < 0x3ffe8000)
 112     {                           /* |x|<0.5 */
 113       if (ix < 0x3fde8000)
 114         {                       /* if |x| < 2**-33 */
 115           if (fabsl (x) < LDBL_MIN)
 116             {
 117               long double force_underflow = x * x;
 118               math_force_eval (force_underflow);
 119             }
 120           if (huge + x > one)
 121             return x;           /* return x with inexact if x!=0 */
 122         }
 123       else
 124         {
 125           t = x * x;
 126           p =
 127             t * (pS0 +
 128                  t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
 129           q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t))));
 130           w = p / q;
 131           return x + x * w;
 132         }
 133     }
 134   /* 1> |x|>= 0.5 */
 135   w = one - fabsl (x);
 136   t = w * 0.5;
 137   p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
 138   q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t))));
 139   s = __ieee754_sqrtl (t);
 140   if (ix >= 0x3ffef999)
 141     {                           /* if |x| > 0.975 */
 142       w = p / q;
 143       t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
 144     }
 145   else
 146     {
 147       GET_LDOUBLE_WORDS (k, i0, i1, s);
 148       i1 = 0;
 149       SET_LDOUBLE_WORDS (w,k,i0,i1);
 150       c = (t - w * w) / (s + w);
 151       r = p / q;
 152       p = 2.0 * s * r - (pio2_lo - 2.0 * c);
 153       q = pio4_hi - 2.0 * w;
 154       t = pio4_hi - (p - q);
 155     }
 156   if ((se & 0x8000) == 0)
 157     return t;
 158   else
 159     return -t;
 160 }
 161 strong_alias (__ieee754_asinl, __asinl_finite)