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