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