sysdeps/ieee754/ldbl-128ibm/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 the
  18   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  *      Between .5 and .625 the approximation is
  40  *              asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
  41  *      For x in [0.625,1]
  42  *              asin(x) = pi/2-2*asin(sqrt((1-x)/2))
  43  *      Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
  44  *      then for x>0.98
  45  *              asin(x) = pi/2 - 2*(s+s*z*R(z))
  46  *                      = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
  47  *      For x<=0.98, let pio4_hi = pio2_hi/2, then
  48  *              f = hi part of s;
  49  *              c = sqrt(z) - f = (z-f*f)/(s+f)         ...f+c=sqrt(z)
  50  *      and
  51  *              asin(x) = pi/2 - 2*(s+s*z*R(z))
  52  *                      = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
  53  *                      = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
  54  *
  55  * Special cases:
  56  *      if x is NaN, return x itself;
  57  *      if |x|>1, return NaN with invalid signal.
  58  *
  59  */
  60
  61
  62 #include "math.h"
  63 #include "math_private.h"
  64 long double sqrtl (long double);
  65
  66 static const long double
  67   one = 1.0L,
  68   huge = 1.0e+300L,
  69   pio2_hi = 1.5707963267948966192313216916397514420986L,
  70   pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
  71   pio4_hi = 7.8539816339744830961566084581987569936977E-1L,
  72
  73         /* coefficient for R(x^2) */
  74
  75   /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
  76      0 <= x <= 0.5
  77      peak relative error 1.9e-35  */
  78   pS0 = -8.358099012470680544198472400254596543711E2L,
  79   pS1 =  3.674973957689619490312782828051860366493E3L,
  80   pS2 = -6.730729094812979665807581609853656623219E3L,
  81   pS3 =  6.643843795209060298375552684423454077633E3L,
  82   pS4 = -3.817341990928606692235481812252049415993E3L,
  83   pS5 =  1.284635388402653715636722822195716476156E3L,
  84   pS6 = -2.410736125231549204856567737329112037867E2L,
  85   pS7 =  2.219191969382402856557594215833622156220E1L,
  86   pS8 = -7.249056260830627156600112195061001036533E-1L,
  87   pS9 =  1.055923570937755300061509030361395604448E-3L,
  88
  89   qS0 = -5.014859407482408326519083440151745519205E3L,
  90   qS1 =  2.430653047950480068881028451580393430537E4L,
  91   qS2 = -4.997904737193653607449250593976069726962E4L,
  92   qS3 =  5.675712336110456923807959930107347511086E4L,
  93   qS4 = -3.881523118339661268482937768522572588022E4L,
  94   qS5 =  1.634202194895541569749717032234510811216E4L,
  95   qS6 = -4.151452662440709301601820849901296953752E3L,
  96   qS7 =  5.956050864057192019085175976175695342168E2L,
  97   qS8 = -4.175375777334867025769346564600396877176E1L,
  98   /* 1.000000000000000000000000000000000000000E0 */
  99
 100   /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
 101      -0.0625 <= x <= 0.0625
 102      peak relative error 3.3e-35  */
 103   rS0 = -5.619049346208901520945464704848780243887E0L,
 104   rS1 =  4.460504162777731472539175700169871920352E1L,
 105   rS2 = -1.317669505315409261479577040530751477488E2L,
 106   rS3 =  1.626532582423661989632442410808596009227E2L,
 107   rS4 = -3.144806644195158614904369445440583873264E1L,
 108   rS5 = -9.806674443470740708765165604769099559553E1L,
 109   rS6 =  5.708468492052010816555762842394927806920E1L,
 110   rS7 =  1.396540499232262112248553357962639431922E1L,
 111   rS8 = -1.126243289311910363001762058295832610344E1L,
 112   rS9 = -4.956179821329901954211277873774472383512E-1L,
 113   rS10 =  3.313227657082367169241333738391762525780E-1L,
 114
 115   sS0 = -4.645814742084009935700221277307007679325E0L,
 116   sS1 =  3.879074822457694323970438316317961918430E1L,
 117   sS2 = -1.221986588013474694623973554726201001066E2L,
 118   sS3 =  1.658821150347718105012079876756201905822E2L,
 119   sS4 = -4.804379630977558197953176474426239748977E1L,
 120   sS5 = -1.004296417397316948114344573811562952793E2L,
 121   sS6 =  7.530281592861320234941101403870010111138E1L,
 122   sS7 =  1.270735595411673647119592092304357226607E1L,
 123   sS8 = -1.815144839646376500705105967064792930282E1L,
 124   sS9 = -7.821597334910963922204235247786840828217E-2L,
 125   /*  1.000000000000000000000000000000000000000E0 */
 126
 127  asinr5625 =  5.9740641664535021430381036628424864397707E-1L;
 128
 129
 130
 131 long double
 132 __ieee754_asinl (long double x)
 133 {
 134   long double t, w, p, q, c, r, s;
 135   int32_t ix, sign, flag;
 136   ieee854_long_double_shape_type u;
 137
 138   flag = 0;
 139   u.value = x;
 140   sign = u.parts32.w0;
 141   ix = sign & 0x7fffffff;
 142   u.parts32.w0 = ix;    /* |x| */
 143   if (ix >= 0x3ff00000) /* |x|>= 1 */
 144     {
 145       if (ix == 0x3ff00000
 146           && (u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3) == 0)
 147         /* asin(1)=+-pi/2 with inexact */
 148         return x * pio2_hi + x * pio2_lo;
 149       return (x - x) / (x - x); /* asin(|x|>1) is NaN */
 150     }
 151   else if (ix < 0x3fe00000) /* |x| < 0.5 */
 152     {
 153       if (ix < 0x3c600000) /* |x| < 2**-57 */
 154         {
 155           if (huge + x > one)
 156             return x;           /* return x with inexact if x!=0 */
 157         }
 158       else
 159         {
 160           t = x * x;
 161           /* Mark to use pS, qS later on.  */
 162           flag = 1;
 163         }
 164     }
 165   else if (ix < 0x3fe40000) /* 0.625 */
 166     {
 167       t = u.value - 0.5625;
 168       p = ((((((((((rS10 * t
 169                     + rS9) * t
 170                    + rS8) * t
 171                   + rS7) * t
 172                  + rS6) * t
 173                 + rS5) * t
 174                + rS4) * t
 175               + rS3) * t
 176              + rS2) * t
 177             + rS1) * t
 178            + rS0) * t;
 179
 180       q = ((((((((( t
 181                     + sS9) * t
 182                   + sS8) * t
 183                  + sS7) * t
 184                 + sS6) * t
 185                + sS5) * t
 186               + sS4) * t
 187              + sS3) * t
 188             + sS2) * t
 189            + sS1) * t
 190         + sS0;
 191       t = asinr5625 + p / q;
 192       if ((sign & 0x80000000) == 0)
 193         return t;
 194       else
 195         return -t;
 196     }
 197   else
 198     {
 199       /* 1 > |x| >= 0.625 */
 200       w = one - u.value;
 201       t = w * 0.5;
 202     }
 203
 204   p = (((((((((pS9 * t
 205                + pS8) * t
 206               + pS7) * t
 207              + pS6) * t
 208             + pS5) * t
 209            + pS4) * t
 210           + pS3) * t
 211          + pS2) * t
 212         + pS1) * t
 213        + pS0) * t;
 214
 215   q = (((((((( t
 216               + qS8) * t
 217              + qS7) * t
 218             + qS6) * t
 219            + qS5) * t
 220           + qS4) * t
 221          + qS3) * t
 222         + qS2) * t
 223        + qS1) * t
 224     + qS0;
 225
 226   if (flag) /* 2^-57 < |x| < 0.5 */
 227     {
 228       w = p / q;
 229       return x + x * w;
 230     }
 231
 232   s = __ieee754_sqrtl (t);
 233   if (ix >= 0x3fef3333) /* |x| > 0.975 */
 234     {
 235       w = p / q;
 236       t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
 237     }
 238   else
 239     {
 240       u.value = s;
 241       u.parts32.w3 = 0;
 242       u.parts32.w2 = 0;
 243       w = u.value;
 244       c = (t - w * w) / (s + w);
 245       r = p / q;
 246       p = 2.0 * s * r - (pio2_lo - 2.0 * c);
 247       q = pio4_hi - 2.0 * w;
 248       t = pio4_hi - (p - q);
 249     }
 250
 251   if ((sign & 0x80000000) == 0)
 252     return t;
 253   else
 254     return -t;
 255 }
 256 strong_alias (__ieee754_asinl, __asinl_finite)