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