Back to index

glibc  2.9
e_asinl.c
Go to the documentation of this file.
00001 /*
00002  * ====================================================
00003  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00004  *
00005  * Developed at SunPro, a Sun Microsystems, Inc. business.
00006  * Permission to use, copy, modify, and distribute this
00007  * software is freely granted, provided that this notice
00008  * is preserved.
00009  * ====================================================
00010  */
00011 
00012 /*
00013   Long double expansions are
00014   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
00015   and are incorporated herein by permission of the author.  The author 
00016   reserves the right to distribute this material elsewhere under different 
00017   copying permissions.  These modifications are distributed here under the 
00018   following terms:
00019 
00020     This library is free software; you can redistribute it and/or
00021     modify it under the terms of the GNU Lesser General Public
00022     License as published by the Free Software Foundation; either
00023     version 2.1 of the License, or (at your option) any later version.
00024 
00025     This library is distributed in the hope that it will be useful,
00026     but WITHOUT ANY WARRANTY; without even the implied warranty of
00027     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00028     Lesser General Public License for more details.
00029 
00030     You should have received a copy of the GNU Lesser General Public
00031     License along with this library; if not, write to the Free Software
00032     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00033 
00034 /* __ieee754_asin(x)
00035  * Method :
00036  *     Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
00037  *     we approximate asin(x) on [0,0.5] by
00038  *            asin(x) = x + x*x^2*R(x^2)
00039  *      Between .5 and .625 the approximation is
00040  *              asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
00041  *     For x in [0.625,1]
00042  *            asin(x) = pi/2-2*asin(sqrt((1-x)/2))
00043  *     Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
00044  *     then for x>0.98
00045  *            asin(x) = pi/2 - 2*(s+s*z*R(z))
00046  *                   = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
00047  *     For x<=0.98, let pio4_hi = pio2_hi/2, then
00048  *            f = hi part of s;
00049  *            c = sqrt(z) - f = (z-f*f)/(s+f)    ...f+c=sqrt(z)
00050  *     and
00051  *            asin(x) = pi/2 - 2*(s+s*z*R(z))
00052  *                   = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
00053  *                   = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
00054  *
00055  * Special cases:
00056  *     if x is NaN, return x itself;
00057  *     if |x|>1, return NaN with invalid signal.
00058  *
00059  */
00060 
00061 
00062 #include "math.h"
00063 #include "math_private.h"
00064 long double sqrtl (long double);
00065 
00066 #ifdef __STDC__
00067 static const long double
00068 #else
00069 static long double
00070 #endif
00071   one = 1.0L,
00072   huge = 1.0e+4932L,
00073   pio2_hi = 1.5707963267948966192313216916397514420986L,
00074   pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
00075   pio4_hi = 7.8539816339744830961566084581987569936977E-1L,
00076 
00077        /* coefficient for R(x^2) */
00078 
00079   /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
00080      0 <= x <= 0.5
00081      peak relative error 1.9e-35  */
00082   pS0 = -8.358099012470680544198472400254596543711E2L,
00083   pS1 =  3.674973957689619490312782828051860366493E3L,
00084   pS2 = -6.730729094812979665807581609853656623219E3L,
00085   pS3 =  6.643843795209060298375552684423454077633E3L,
00086   pS4 = -3.817341990928606692235481812252049415993E3L,
00087   pS5 =  1.284635388402653715636722822195716476156E3L,
00088   pS6 = -2.410736125231549204856567737329112037867E2L,
00089   pS7 =  2.219191969382402856557594215833622156220E1L,
00090   pS8 = -7.249056260830627156600112195061001036533E-1L,
00091   pS9 =  1.055923570937755300061509030361395604448E-3L,
00092 
00093   qS0 = -5.014859407482408326519083440151745519205E3L,
00094   qS1 =  2.430653047950480068881028451580393430537E4L,
00095   qS2 = -4.997904737193653607449250593976069726962E4L,
00096   qS3 =  5.675712336110456923807959930107347511086E4L,
00097   qS4 = -3.881523118339661268482937768522572588022E4L,
00098   qS5 =  1.634202194895541569749717032234510811216E4L,
00099   qS6 = -4.151452662440709301601820849901296953752E3L,
00100   qS7 =  5.956050864057192019085175976175695342168E2L,
00101   qS8 = -4.175375777334867025769346564600396877176E1L,
00102   /* 1.000000000000000000000000000000000000000E0 */
00103 
00104   /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
00105      -0.0625 <= x <= 0.0625
00106      peak relative error 3.3e-35  */
00107   rS0 = -5.619049346208901520945464704848780243887E0L,
00108   rS1 =  4.460504162777731472539175700169871920352E1L,
00109   rS2 = -1.317669505315409261479577040530751477488E2L,
00110   rS3 =  1.626532582423661989632442410808596009227E2L,
00111   rS4 = -3.144806644195158614904369445440583873264E1L,
00112   rS5 = -9.806674443470740708765165604769099559553E1L,
00113   rS6 =  5.708468492052010816555762842394927806920E1L,
00114   rS7 =  1.396540499232262112248553357962639431922E1L,
00115   rS8 = -1.126243289311910363001762058295832610344E1L,
00116   rS9 = -4.956179821329901954211277873774472383512E-1L,
00117   rS10 =  3.313227657082367169241333738391762525780E-1L,
00118 
00119   sS0 = -4.645814742084009935700221277307007679325E0L,
00120   sS1 =  3.879074822457694323970438316317961918430E1L,
00121   sS2 = -1.221986588013474694623973554726201001066E2L,
00122   sS3 =  1.658821150347718105012079876756201905822E2L,
00123   sS4 = -4.804379630977558197953176474426239748977E1L,
00124   sS5 = -1.004296417397316948114344573811562952793E2L,
00125   sS6 =  7.530281592861320234941101403870010111138E1L,
00126   sS7 =  1.270735595411673647119592092304357226607E1L,
00127   sS8 = -1.815144839646376500705105967064792930282E1L,
00128   sS9 = -7.821597334910963922204235247786840828217E-2L,
00129   /*  1.000000000000000000000000000000000000000E0 */
00130 
00131  asinr5625 =  5.9740641664535021430381036628424864397707E-1L;
00132 
00133 
00134 
00135 #ifdef __STDC__
00136 long double
00137 __ieee754_asinl (long double x)
00138 #else
00139 double
00140 __ieee754_asinl (x)
00141      long double x;
00142 #endif
00143 {
00144   long double t, w, p, q, c, r, s;
00145   int32_t ix, sign, flag;
00146   ieee854_long_double_shape_type u;
00147 
00148   flag = 0;
00149   u.value = x;
00150   sign = u.parts32.w0;
00151   ix = sign & 0x7fffffff;
00152   u.parts32.w0 = ix;    /* |x| */
00153   if (ix >= 0x3fff0000)     /* |x|>= 1 */
00154     {
00155       if (ix == 0x3fff0000
00156          && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
00157        /* asin(1)=+-pi/2 with inexact */
00158        return x * pio2_hi + x * pio2_lo;
00159       return (x - x) / (x - x);    /* asin(|x|>1) is NaN */
00160     }
00161   else if (ix < 0x3ffe0000) /* |x| < 0.5 */
00162     {
00163       if (ix < 0x3fc60000) /* |x| < 2**-57 */
00164        {
00165          if (huge + x > one)
00166            return x;        /* return x with inexact if x!=0 */
00167        }
00168       else
00169        {
00170          t = x * x;
00171          /* Mark to use pS, qS later on.  */
00172          flag = 1;
00173        }
00174     }
00175   else if (ix < 0x3ffe4000) /* 0.625 */
00176     {
00177       t = u.value - 0.5625;
00178       p = ((((((((((rS10 * t
00179                   + rS9) * t
00180                  + rS8) * t
00181                 + rS7) * t
00182                + rS6) * t
00183               + rS5) * t
00184               + rS4) * t
00185              + rS3) * t
00186             + rS2) * t
00187            + rS1) * t
00188           + rS0) * t;
00189 
00190       q = ((((((((( t
00191                   + sS9) * t
00192                 + sS8) * t
00193                + sS7) * t
00194               + sS6) * t
00195               + sS5) * t
00196              + sS4) * t
00197             + sS3) * t
00198            + sS2) * t
00199           + sS1) * t
00200        + sS0;
00201       t = asinr5625 + p / q;
00202       if ((sign & 0x80000000) == 0)
00203        return t;
00204       else
00205        return -t;
00206     }
00207   else
00208     {
00209       /* 1 > |x| >= 0.625 */
00210       w = one - u.value;
00211       t = w * 0.5;
00212     }
00213 
00214   p = (((((((((pS9 * t
00215               + pS8) * t
00216              + pS7) * t
00217             + pS6) * t
00218            + pS5) * t
00219           + pS4) * t
00220          + pS3) * t
00221         + pS2) * t
00222        + pS1) * t
00223        + pS0) * t;
00224 
00225   q = (((((((( t
00226              + qS8) * t
00227             + qS7) * t
00228            + qS6) * t
00229           + qS5) * t
00230          + qS4) * t
00231         + qS3) * t
00232        + qS2) * t
00233        + qS1) * t
00234     + qS0;
00235 
00236   if (flag) /* 2^-57 < |x| < 0.5 */
00237     {
00238       w = p / q;
00239       return x + x * w;
00240     }
00241 
00242   s = __ieee754_sqrtl (t);
00243   if (ix >= 0x3ffef333) /* |x| > 0.975 */
00244     {
00245       w = p / q;
00246       t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
00247     }
00248   else
00249     {
00250       u.value = s;
00251       u.parts32.w3 = 0;
00252       u.parts32.w2 = 0;
00253       w = u.value;
00254       c = (t - w * w) / (s + w);
00255       r = p / q;
00256       p = 2.0 * s * r - (pio2_lo - 2.0 * c);
00257       q = pio4_hi - 2.0 * w;
00258       t = pio4_hi - (p - q);
00259     }
00260 
00261   if ((sign & 0x80000000) == 0)
00262     return t;
00263   else
00264     return -t;
00265 }