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     <https://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 #include <math-underflow.h>
  65
  66 static const long double
  67   one = 1.0L,
  68   huge = 1.0e+4932L,
  69   pio2_hi = 0x1.921fb54442d1846ap+0L, /* pi/2 rounded to nearest to 64
  70                                          bits.  */
  71   pio2_lo = -0x7.6733ae8fe47c65d8p-68L, /* pi/2 - pio2_hi rounded to
  72                                            nearest to 64 bits.  */
  73   pio4_hi = 0xc.90fdaa22168c235p-4L, /* pi/4 rounded to nearest to 64
  74                                         bits.  */
  75
  76         /* coefficient for R(x^2) */
  77
  78   /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
  79      0 <= x <= 0.5
  80      peak relative error 1.9e-21  */
  81   pS0 =  -1.008714657938491626019651170502036851607E1L,
  82   pS1 =   2.331460313214179572063441834101394865259E1L,
  83   pS2 =  -1.863169762159016144159202387315381830227E1L,
  84   pS3 =   5.930399351579141771077475766877674661747E0L,
  85   pS4 =  -6.121291917696920296944056882932695185001E-1L,
  86   pS5 =   3.776934006243367487161248678019350338383E-3L,
  87
  88   qS0 =  -6.052287947630949712886794360635592886517E1L,
  89   qS1 =   1.671229145571899593737596543114258558503E2L,
  90   qS2 =  -1.707840117062586426144397688315411324388E2L,
  91   qS3 =   7.870295154902110425886636075950077640623E1L,
  92   qS4 =  -1.568433562487314651121702982333303458814E1L;
  93     /* 1.000000000000000000000000000000000000000E0 */
  94
  95 long double
  96 __ieee754_asinl (long double x)
  97 {
  98   long double t, w, p, q, c, r, s;
  99   int32_t ix;
 100   uint32_t se, i0, i1, k;
 101
 102   GET_LDOUBLE_WORDS (se, i0, i1, x);
 103   ix = se & 0x7fff;
 104   ix = (ix << 16) | (i0 >> 16);
 105   if (ix >= 0x3fff8000)
 106     {                           /* |x|>= 1 */
 107       if (ix == 0x3fff8000 && ((i0 - 0x80000000) | i1) == 0)
 108         /* asin(1)=+-pi/2 with inexact */
 109         return x * pio2_hi + x * pio2_lo;
 110       return (x - x) / (x - x); /* asin(|x|>1) is NaN */
 111     }
 112   else if (ix < 0x3ffe8000)
 113     {                           /* |x|<0.5 */
 114       if (ix < 0x3fde8000)
 115         {                       /* if |x| < 2**-33 */
 116           math_check_force_underflow (x);
 117           if (huge + x > one)
 118             return x;           /* return x with inexact if x!=0 */
 119         }
 120       else
 121         {
 122           t = x * x;
 123           p =
 124             t * (pS0 +
 125                  t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
 126           q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t))));
 127           w = p / q;
 128           return x + x * w;
 129         }
 130     }
 131   /* 1> |x|>= 0.5 */
 132   w = one - fabsl (x);
 133   t = w * 0.5;
 134   p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
 135   q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t))));
 136   s = sqrtl (t);
 137   if (ix >= 0x3ffef999)
 138     {                           /* if |x| > 0.975 */
 139       w = p / q;
 140       t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
 141     }
 142   else
 143     {
 144       GET_LDOUBLE_WORDS (k, i0, i1, s);
 145       i1 = 0;
 146       SET_LDOUBLE_WORDS (w,k,i0,i1);
 147       c = (t - w * w) / (s + w);
 148       r = p / q;
 149       p = 2.0 * s * r - (pio2_lo - 2.0 * c);
 150       q = pio4_hi - 2.0 * w;
 151       t = pio4_hi - (p - q);
 152     }
 153   if ((se & 0x8000) == 0)
 154     return t;
 155   else
 156     return -t;
 157 }
 158 strong_alias (__ieee754_asinl, __asinl_finite)