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