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