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, write to the Free Software
  32     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
  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 #ifdef __STDC__
  67 static const long double
  68 #else
  69 static long double
  70 #endif
  71   one = 1.0L,
  72   huge = 1.0e+4932L,
  73   pio2_hi = 1.5707963267948966192313216916397514420986L,
  74   pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
  75   pio4_hi = 7.8539816339744830961566084581987569936977E-1L,
  76
  77         /* coefficient for R(x^2) */
  78
  79   /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
  80      0 <= x <= 0.5
  81      peak relative error 1.9e-35  */
  82   pS0 = -8.358099012470680544198472400254596543711E2L,
  83   pS1 =  3.674973957689619490312782828051860366493E3L,
  84   pS2 = -6.730729094812979665807581609853656623219E3L,
  85   pS3 =  6.643843795209060298375552684423454077633E3L,
  86   pS4 = -3.817341990928606692235481812252049415993E3L,
  87   pS5 =  1.284635388402653715636722822195716476156E3L,
  88   pS6 = -2.410736125231549204856567737329112037867E2L,
  89   pS7 =  2.219191969382402856557594215833622156220E1L,
  90   pS8 = -7.249056260830627156600112195061001036533E-1L,
  91   pS9 =  1.055923570937755300061509030361395604448E-3L,
  92
  93   qS0 = -5.014859407482408326519083440151745519205E3L,
  94   qS1 =  2.430653047950480068881028451580393430537E4L,
  95   qS2 = -4.997904737193653607449250593976069726962E4L,
  96   qS3 =  5.675712336110456923807959930107347511086E4L,
  97   qS4 = -3.881523118339661268482937768522572588022E4L,
  98   qS5 =  1.634202194895541569749717032234510811216E4L,
  99   qS6 = -4.151452662440709301601820849901296953752E3L,
 100   qS7 =  5.956050864057192019085175976175695342168E2L,
 101   qS8 = -4.175375777334867025769346564600396877176E1L,
 102   /* 1.000000000000000000000000000000000000000E0 */
 103
 104   /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
 105      -0.0625 <= x <= 0.0625
 106      peak relative error 3.3e-35  */
 107   rS0 = -5.619049346208901520945464704848780243887E0L,
 108   rS1 =  4.460504162777731472539175700169871920352E1L,
 109   rS2 = -1.317669505315409261479577040530751477488E2L,
 110   rS3 =  1.626532582423661989632442410808596009227E2L,
 111   rS4 = -3.144806644195158614904369445440583873264E1L,
 112   rS5 = -9.806674443470740708765165604769099559553E1L,
 113   rS6 =  5.708468492052010816555762842394927806920E1L,
 114   rS7 =  1.396540499232262112248553357962639431922E1L,
 115   rS8 = -1.126243289311910363001762058295832610344E1L,
 116   rS9 = -4.956179821329901954211277873774472383512E-1L,
 117   rS10 =  3.313227657082367169241333738391762525780E-1L,
 118
 119   sS0 = -4.645814742084009935700221277307007679325E0L,
 120   sS1 =  3.879074822457694323970438316317961918430E1L,
 121   sS2 = -1.221986588013474694623973554726201001066E2L,
 122   sS3 =  1.658821150347718105012079876756201905822E2L,
 123   sS4 = -4.804379630977558197953176474426239748977E1L,
 124   sS5 = -1.004296417397316948114344573811562952793E2L,
 125   sS6 =  7.530281592861320234941101403870010111138E1L,
 126   sS7 =  1.270735595411673647119592092304357226607E1L,
 127   sS8 = -1.815144839646376500705105967064792930282E1L,
 128   sS9 = -7.821597334910963922204235247786840828217E-2L,
 129   /*  1.000000000000000000000000000000000000000E0 */
 130
 131  asinr5625 =  5.9740641664535021430381036628424864397707E-1L;
 132
 133
 134
 135 #ifdef __STDC__
 136 long double
 137 __ieee754_asinl (long double x)
 138 #else
 139 double
 140 __ieee754_asinl (x)
 141      long double x;
 142 #endif
 143 {
 144   long double t, w, p, q, c, r, s;
 145   int32_t ix, sign, flag;
 146   ieee854_long_double_shape_type u;
 147
 148   flag = 0;
 149   u.value = x;
 150   sign = u.parts32.w0;
 151   ix = sign & 0x7fffffff;
 152   u.parts32.w0 = ix;    /* |x| */
 153   if (ix >= 0x3fff0000) /* |x|>= 1 */
 154     {
 155       if (ix == 0x3fff0000
 156           && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
 157         /* asin(1)=+-pi/2 with inexact */
 158         return x * pio2_hi + x * pio2_lo;
 159       return (x - x) / (x - x); /* asin(|x|>1) is NaN */
 160     }
 161   else if (ix < 0x3ffe0000) /* |x| < 0.5 */
 162     {
 163       if (ix < 0x3fc60000) /* |x| < 2**-57 */
 164         {
 165           if (huge + x > one)
 166             return x;           /* return x with inexact if x!=0 */
 167         }
 168       else
 169         {
 170           t = x * x;
 171           /* Mark to use pS, qS later on.  */
 172           flag = 1;
 173         }
 174     }
 175   else if (ix < 0x3ffe4000) /* 0.625 */
 176     {
 177       t = u.value - 0.5625;
 178       p = ((((((((((rS10 * t
 179                     + rS9) * t
 180                    + rS8) * t
 181                   + rS7) * t
 182                  + rS6) * t
 183                 + rS5) * t
 184                + rS4) * t
 185               + rS3) * t
 186              + rS2) * t
 187             + rS1) * t
 188            + rS0) * t;
 189
 190       q = ((((((((( t
 191                     + sS9) * t
 192                   + sS8) * t
 193                  + sS7) * t
 194                 + sS6) * t
 195                + sS5) * t
 196               + sS4) * t
 197              + sS3) * t
 198             + sS2) * t
 199            + sS1) * t
 200         + sS0;
 201       t = asinr5625 + p / q;
 202       if ((sign & 0x80000000) == 0)
 203         return t;
 204       else
 205         return -t;
 206     }
 207   else
 208     {
 209       /* 1 > |x| >= 0.625 */
 210       w = one - u.value;
 211       t = w * 0.5;
 212     }
 213
 214   p = (((((((((pS9 * t
 215                + pS8) * t
 216               + pS7) * t
 217              + pS6) * t
 218             + pS5) * t
 219            + pS4) * t
 220           + pS3) * t
 221          + pS2) * t
 222         + pS1) * t
 223        + pS0) * t;
 224
 225   q = (((((((( t
 226               + qS8) * t
 227              + qS7) * t
 228             + qS6) * t
 229            + qS5) * t
 230           + qS4) * t
 231          + qS3) * t
 232         + qS2) * t
 233        + qS1) * t
 234     + qS0;
 235
 236   if (flag) /* 2^-57 < |x| < 0.5 */
 237     {
 238       w = p / q;
 239       return x + x * w;
 240     }
 241
 242   s = __ieee754_sqrtl (t);
 243   if (ix >= 0x3ffef333) /* |x| > 0.975 */
 244     {
 245       w = p / q;
 246       t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
 247     }
 248   else
 249     {
 250       u.value = s;
 251       u.parts32.w3 = 0;
 252       u.parts32.w2 = 0;
 253       w = u.value;
 254       c = (t - w * w) / (s + w);
 255       r = p / q;
 256       p = 2.0 * s * r - (pio2_lo - 2.0 * c);
 257       q = pio4_hi - 2.0 * w;
 258       t = pio4_hi - (p - q);
 259     }
 260
 261   if ((sign & 0x80000000) == 0)
 262     return t;
 263   else
 264     return -t;
 265 }