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           if (huge + x > one)
 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 }