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 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 <float.h>
  63 #include <math.h>
  64 #include <math_private.h>
  65
  66 static const _Float128
  67   one = 1,
  68   huge = L(1.0e+4932),
  69   pio2_hi = L(1.5707963267948966192313216916397514420986),
  70   pio2_lo = L(4.3359050650618905123985220130216759843812E-35),
  71   pio4_hi = L(7.8539816339744830961566084581987569936977E-1),
  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 = L(-8.358099012470680544198472400254596543711E2),
  79   pS1 =  L(3.674973957689619490312782828051860366493E3),
  80   pS2 = L(-6.730729094812979665807581609853656623219E3),
  81   pS3 =  L(6.643843795209060298375552684423454077633E3),
  82   pS4 = L(-3.817341990928606692235481812252049415993E3),
  83   pS5 =  L(1.284635388402653715636722822195716476156E3),
  84   pS6 = L(-2.410736125231549204856567737329112037867E2),
  85   pS7 =  L(2.219191969382402856557594215833622156220E1),
  86   pS8 = L(-7.249056260830627156600112195061001036533E-1),
  87   pS9 =  L(1.055923570937755300061509030361395604448E-3),
  88
  89   qS0 = L(-5.014859407482408326519083440151745519205E3),
  90   qS1 =  L(2.430653047950480068881028451580393430537E4),
  91   qS2 = L(-4.997904737193653607449250593976069726962E4),
  92   qS3 =  L(5.675712336110456923807959930107347511086E4),
  93   qS4 = L(-3.881523118339661268482937768522572588022E4),
  94   qS5 =  L(1.634202194895541569749717032234510811216E4),
  95   qS6 = L(-4.151452662440709301601820849901296953752E3),
  96   qS7 =  L(5.956050864057192019085175976175695342168E2),
  97   qS8 = L(-4.175375777334867025769346564600396877176E1),
  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 = L(-5.619049346208901520945464704848780243887E0),
 104   rS1 =  L(4.460504162777731472539175700169871920352E1),
 105   rS2 = L(-1.317669505315409261479577040530751477488E2),
 106   rS3 =  L(1.626532582423661989632442410808596009227E2),
 107   rS4 = L(-3.144806644195158614904369445440583873264E1),
 108   rS5 = L(-9.806674443470740708765165604769099559553E1),
 109   rS6 =  L(5.708468492052010816555762842394927806920E1),
 110   rS7 =  L(1.396540499232262112248553357962639431922E1),
 111   rS8 = L(-1.126243289311910363001762058295832610344E1),
 112   rS9 = L(-4.956179821329901954211277873774472383512E-1),
 113   rS10 =  L(3.313227657082367169241333738391762525780E-1),
 114
 115   sS0 = L(-4.645814742084009935700221277307007679325E0),
 116   sS1 =  L(3.879074822457694323970438316317961918430E1),
 117   sS2 = L(-1.221986588013474694623973554726201001066E2),
 118   sS3 =  L(1.658821150347718105012079876756201905822E2),
 119   sS4 = L(-4.804379630977558197953176474426239748977E1),
 120   sS5 = L(-1.004296417397316948114344573811562952793E2),
 121   sS6 =  L(7.530281592861320234941101403870010111138E1),
 122   sS7 =  L(1.270735595411673647119592092304357226607E1),
 123   sS8 = L(-1.815144839646376500705105967064792930282E1),
 124   sS9 = L(-7.821597334910963922204235247786840828217E-2),
 125   /*  1.000000000000000000000000000000000000000E0 */
 126
 127  asinr5625 =  L(5.9740641664535021430381036628424864397707E-1);
 128
 129
 130
 131 _Float128
 132 __ieee754_asinl (_Float128 x)
 133 {
 134   _Float128 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 >= 0x3fff0000) /* |x|>= 1 */
 144     {
 145       if (ix == 0x3fff0000
 146           && (u.parts32.w1 | u.parts32.w2 | 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 < 0x3ffe0000) /* |x| < 0.5 */
 152     {
 153       if (ix < 0x3fc60000) /* |x| < 2**-57 */
 154         {
 155           math_check_force_underflow (x);
 156           _Float128 force_inexact = huge + x;
 157           math_force_eval (force_inexact);
 158           return x;             /* return x with inexact if x!=0 */
 159         }
 160       else
 161         {
 162           t = x * x;
 163           /* Mark to use pS, qS later on.  */
 164           flag = 1;
 165         }
 166     }
 167   else if (ix < 0x3ffe4000) /* 0.625 */
 168     {
 169       t = u.value - 0.5625;
 170       p = ((((((((((rS10 * t
 171                     + rS9) * t
 172                    + rS8) * t
 173                   + rS7) * t
 174                  + rS6) * t
 175                 + rS5) * t
 176                + rS4) * t
 177               + rS3) * t
 178              + rS2) * t
 179             + rS1) * t
 180            + rS0) * t;
 181
 182       q = ((((((((( t
 183                     + sS9) * t
 184                   + sS8) * t
 185                  + sS7) * t
 186                 + sS6) * t
 187                + sS5) * t
 188               + sS4) * t
 189              + sS3) * t
 190             + sS2) * t
 191            + sS1) * t
 192         + sS0;
 193       t = asinr5625 + p / q;
 194       if ((sign & 0x80000000) == 0)
 195         return t;
 196       else
 197         return -t;
 198     }
 199   else
 200     {
 201       /* 1 > |x| >= 0.625 */
 202       w = one - u.value;
 203       t = w * 0.5;
 204     }
 205
 206   p = (((((((((pS9 * t
 207                + pS8) * t
 208               + pS7) * t
 209              + pS6) * t
 210             + pS5) * t
 211            + pS4) * t
 212           + pS3) * t
 213          + pS2) * t
 214         + pS1) * t
 215        + pS0) * t;
 216
 217   q = (((((((( t
 218               + qS8) * t
 219              + qS7) * t
 220             + qS6) * t
 221            + qS5) * t
 222           + qS4) * t
 223          + qS3) * t
 224         + qS2) * t
 225        + qS1) * t
 226     + qS0;
 227
 228   if (flag) /* 2^-57 < |x| < 0.5 */
 229     {
 230       w = p / q;
 231       return x + x * w;
 232     }
 233
 234   s = __ieee754_sqrtl (t);
 235   if (ix >= 0x3ffef333) /* |x| > 0.975 */
 236     {
 237       w = p / q;
 238       t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
 239     }
 240   else
 241     {
 242       u.value = s;
 243       u.parts32.w3 = 0;
 244       u.parts32.w2 = 0;
 245       w = u.value;
 246       c = (t - w * w) / (s + w);
 247       r = p / q;
 248       p = 2.0 * s * r - (pio2_lo - 2.0 * c);
 249       q = pio4_hi - 2.0 * w;
 250       t = pio4_hi - (p - q);
 251     }
 252
 253   if ((sign & 0x80000000) == 0)
 254     return t;
 255   else
 256     return -t;
 257 }
 258 strong_alias (__ieee754_asinl, __asinl_finite)