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