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 a, t, w, p, q, c, r, s;
 135   int flag;
 136
 137   if (__glibc_unlikely (__isnanl (x)))
 138     return x + x;
 139   flag = 0;
 140   a = __builtin_fabsl (x);
 141   if (a == 1.0L)        /* |x|>= 1 */
 142     return x * pio2_hi + x * pio2_lo;   /* asin(1)=+-pi/2 with inexact */
 143   else if (a >= 1.0L)
 144     return (x - x) / (x - x);   /* asin(|x|>1) is NaN */
 145   else if (a < 0.5L)
 146     {
 147       if (a < 6.938893903907228e-18L) /* |x| < 2**-57 */
 148         {
 149           if (huge + x > one)
 150             return x;           /* return x with inexact if x!=0 */
 151         }
 152       else
 153         {
 154           t = x * x;
 155           /* Mark to use pS, qS later on.  */
 156           flag = 1;
 157         }
 158     }
 159   else if (a < 0.625L)
 160     {
 161       t = a - 0.5625;
 162       p = ((((((((((rS10 * t
 163                     + rS9) * t
 164                    + rS8) * t
 165                   + rS7) * t
 166                  + rS6) * t
 167                 + rS5) * t
 168                + rS4) * t
 169               + rS3) * t
 170              + rS2) * t
 171             + rS1) * t
 172            + rS0) * t;
 173
 174       q = ((((((((( t
 175                     + sS9) * t
 176                   + sS8) * t
 177                  + sS7) * t
 178                 + sS6) * t
 179                + sS5) * t
 180               + sS4) * t
 181              + sS3) * t
 182             + sS2) * t
 183            + sS1) * t
 184         + sS0;
 185       t = asinr5625 + p / q;
 186       if (x > 0.0L)
 187         return t;
 188       else
 189         return -t;
 190     }
 191   else
 192     {
 193       /* 1 > |x| >= 0.625 */
 194       w = one - a;
 195       t = w * 0.5;
 196     }
 197
 198   p = (((((((((pS9 * t
 199                + pS8) * t
 200               + pS7) * t
 201              + pS6) * t
 202             + pS5) * t
 203            + pS4) * t
 204           + pS3) * t
 205          + pS2) * t
 206         + pS1) * t
 207        + pS0) * t;
 208
 209   q = (((((((( t
 210               + qS8) * t
 211              + qS7) * t
 212             + qS6) * t
 213            + qS5) * t
 214           + qS4) * t
 215          + qS3) * t
 216         + qS2) * t
 217        + qS1) * t
 218     + qS0;
 219
 220   if (flag) /* 2^-57 < |x| < 0.5 */
 221     {
 222       w = p / q;
 223       return x + x * w;
 224     }
 225
 226   s = __ieee754_sqrtl (t);
 227   if (a > 0.975L)
 228     {
 229       w = p / q;
 230       t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
 231     }
 232   else
 233     {
 234       w = ldbl_high (s);
 235       c = (t - w * w) / (s + w);
 236       r = p / q;
 237       p = 2.0 * s * r - (pio2_lo - 2.0 * c);
 238       q = pio4_hi - 2.0 * w;
 239       t = pio4_hi - (p - q);
 240     }
 241
 242   if (x > 0.0L)
 243     return t;
 244   else
 245     return -t;
 246 }
 247 strong_alias (__ieee754_asinl, __asinl_finite)