Back to index

glibc  2.9
e_acosl.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 
00018    the 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_acosl(x)
00035  * Method :
00036  *      acos(x)  = pi/2 - asin(x)
00037  *      acos(-x) = pi/2 + asin(x)
00038  * For |x| <= 0.375
00039  *      acos(x) = pi/2 - asin(x)
00040  * Between .375 and .5 the approximation is
00041  *      acos(0.4375 + x) = acos(0.4375) + x P(x) / Q(x)
00042  * Between .5 and .625 the approximation is
00043  *      acos(0.5625 + x) = acos(0.5625) + x rS(x) / sS(x)
00044  * For x > 0.625,
00045  *      acos(x) = 2 asin(sqrt((1-x)/2))
00046  *      computed with an extended precision square root in the leading term.
00047  * For x < -0.625
00048  *      acos(x) = pi - 2 asin(sqrt((1-|x|)/2))
00049  *
00050  * Special cases:
00051  *      if x is NaN, return x itself;
00052  *      if |x|>1, return NaN with invalid signal.
00053  *
00054  * Functions needed: __ieee754_sqrtl.
00055  */
00056 
00057 #include "math.h"
00058 #include "math_private.h"
00059 
00060 #ifdef __STDC__
00061 static const long double
00062 #else
00063 static long double
00064 #endif
00065   one = 1.0L,
00066   pio2_hi = 1.5707963267948966192313216916397514420986L,
00067   pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
00068 
00069   /* acos(0.5625 + x) = acos(0.5625) + x rS(x) / sS(x)
00070      -0.0625 <= x <= 0.0625
00071      peak relative error 3.3e-35  */
00072 
00073   rS0 =  5.619049346208901520945464704848780243887E0L,
00074   rS1 = -4.460504162777731472539175700169871920352E1L,
00075   rS2 =  1.317669505315409261479577040530751477488E2L,
00076   rS3 = -1.626532582423661989632442410808596009227E2L,
00077   rS4 =  3.144806644195158614904369445440583873264E1L,
00078   rS5 =  9.806674443470740708765165604769099559553E1L,
00079   rS6 = -5.708468492052010816555762842394927806920E1L,
00080   rS7 = -1.396540499232262112248553357962639431922E1L,
00081   rS8 =  1.126243289311910363001762058295832610344E1L,
00082   rS9 =  4.956179821329901954211277873774472383512E-1L,
00083   rS10 = -3.313227657082367169241333738391762525780E-1L,
00084 
00085   sS0 = -4.645814742084009935700221277307007679325E0L,
00086   sS1 =  3.879074822457694323970438316317961918430E1L,
00087   sS2 = -1.221986588013474694623973554726201001066E2L,
00088   sS3 =  1.658821150347718105012079876756201905822E2L,
00089   sS4 = -4.804379630977558197953176474426239748977E1L,
00090   sS5 = -1.004296417397316948114344573811562952793E2L,
00091   sS6 =  7.530281592861320234941101403870010111138E1L,
00092   sS7 =  1.270735595411673647119592092304357226607E1L,
00093   sS8 = -1.815144839646376500705105967064792930282E1L,
00094   sS9 = -7.821597334910963922204235247786840828217E-2L,
00095   /* 1.000000000000000000000000000000000000000E0 */
00096 
00097   acosr5625 = 9.7338991014954640492751132535550279812151E-1L,
00098   pimacosr5625 = 2.1682027434402468335351320579240000860757E0L,
00099 
00100   /* acos(0.4375 + x) = acos(0.4375) + x rS(x) / sS(x)
00101      -0.0625 <= x <= 0.0625
00102      peak relative error 2.1e-35  */
00103 
00104   P0 =  2.177690192235413635229046633751390484892E0L,
00105   P1 = -2.848698225706605746657192566166142909573E1L,
00106   P2 =  1.040076477655245590871244795403659880304E2L,
00107   P3 = -1.400087608918906358323551402881238180553E2L,
00108   P4 =  2.221047917671449176051896400503615543757E1L,
00109   P5 =  9.643714856395587663736110523917499638702E1L,
00110   P6 = -5.158406639829833829027457284942389079196E1L,
00111   P7 = -1.578651828337585944715290382181219741813E1L,
00112   P8 =  1.093632715903802870546857764647931045906E1L,
00113   P9 =  5.448925479898460003048760932274085300103E-1L,
00114   P10 = -3.315886001095605268470690485170092986337E-1L,
00115   Q0 = -1.958219113487162405143608843774587557016E0L,
00116   Q1 =  2.614577866876185080678907676023269360520E1L,
00117   Q2 = -9.990858606464150981009763389881793660938E1L,
00118   Q3 =  1.443958741356995763628660823395334281596E2L,
00119   Q4 = -3.206441012484232867657763518369723873129E1L,
00120   Q5 = -1.048560885341833443564920145642588991492E2L,
00121   Q6 =  6.745883931909770880159915641984874746358E1L,
00122   Q7 =  1.806809656342804436118449982647641392951E1L,
00123   Q8 = -1.770150690652438294290020775359580915464E1L,
00124   Q9 = -5.659156469628629327045433069052560211164E-1L,
00125   /* 1.000000000000000000000000000000000000000E0 */
00126 
00127   acosr4375 = 1.1179797320499710475919903296900511518755E0L,
00128   pimacosr4375 = 2.0236129215398221908706530535894517323217E0L,
00129 
00130   /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
00131      0 <= x <= 0.5
00132      peak relative error 1.9e-35  */
00133   pS0 = -8.358099012470680544198472400254596543711E2L,
00134   pS1 =  3.674973957689619490312782828051860366493E3L,
00135   pS2 = -6.730729094812979665807581609853656623219E3L,
00136   pS3 =  6.643843795209060298375552684423454077633E3L,
00137   pS4 = -3.817341990928606692235481812252049415993E3L,
00138   pS5 =  1.284635388402653715636722822195716476156E3L,
00139   pS6 = -2.410736125231549204856567737329112037867E2L,
00140   pS7 =  2.219191969382402856557594215833622156220E1L,
00141   pS8 = -7.249056260830627156600112195061001036533E-1L,
00142   pS9 =  1.055923570937755300061509030361395604448E-3L,
00143 
00144   qS0 = -5.014859407482408326519083440151745519205E3L,
00145   qS1 =  2.430653047950480068881028451580393430537E4L,
00146   qS2 = -4.997904737193653607449250593976069726962E4L,
00147   qS3 =  5.675712336110456923807959930107347511086E4L,
00148   qS4 = -3.881523118339661268482937768522572588022E4L,
00149   qS5 =  1.634202194895541569749717032234510811216E4L,
00150   qS6 = -4.151452662440709301601820849901296953752E3L,
00151   qS7 =  5.956050864057192019085175976175695342168E2L,
00152   qS8 = -4.175375777334867025769346564600396877176E1L;
00153   /* 1.000000000000000000000000000000000000000E0 */
00154 
00155 #ifdef __STDC__
00156 long double
00157 __ieee754_acosl (long double x)
00158 #else
00159 long double
00160 __ieee754_acosl (x)
00161      long double x;
00162 #endif
00163 {
00164   long double z, r, w, p, q, s, t, f2;
00165   int32_t ix, sign;
00166   ieee854_long_double_shape_type u;
00167 
00168   u.value = x;
00169   sign = u.parts32.w0;
00170   ix = sign & 0x7fffffff;
00171   u.parts32.w0 = ix;        /* |x| */
00172   if (ix >= 0x3fff0000)            /* |x| >= 1 */
00173     {
00174       if (ix == 0x3fff0000
00175          && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
00176        {                    /* |x| == 1 */
00177          if ((sign & 0x80000000) == 0)
00178            return 0.0;             /* acos(1) = 0  */
00179          else
00180            return (2.0 * pio2_hi) + (2.0 * pio2_lo);    /* acos(-1)= pi */
00181        }
00182       return (x - x) / (x - x);    /* acos(|x| > 1) is NaN */
00183     }
00184   else if (ix < 0x3ffe0000) /* |x| < 0.5 */
00185     {
00186       if (ix < 0x3fc60000)  /* |x| < 2**-57 */
00187        return pio2_hi + pio2_lo;
00188       if (ix < 0x3ffde000)  /* |x| < .4375 */
00189        {
00190          /* Arcsine of x.  */
00191          z = x * x;
00192          p = (((((((((pS9 * z
00193                      + pS8) * z
00194                     + pS7) * z
00195                    + pS6) * z
00196                   + pS5) * z
00197                  + pS4) * z
00198                 + pS3) * z
00199                + pS2) * z
00200               + pS1) * z
00201               + pS0) * z;
00202          q = (((((((( z
00203                      + qS8) * z
00204                    + qS7) * z
00205                   + qS6) * z
00206                  + qS5) * z
00207                 + qS4) * z
00208                + qS3) * z
00209               + qS2) * z
00210               + qS1) * z
00211            + qS0;
00212          r = x + x * p / q;
00213          z = pio2_hi - (r - pio2_lo);
00214          return z;
00215        }
00216       /* .4375 <= |x| < .5 */
00217       t = u.value - 0.4375L;
00218       p = ((((((((((P10 * t
00219                   + P9) * t
00220                  + P8) * t
00221                 + P7) * t
00222                + P6) * t
00223               + P5) * t
00224               + P4) * t
00225              + P3) * t
00226             + P2) * t
00227            + P1) * t
00228           + P0) * t;
00229 
00230       q = (((((((((t
00231                  + Q9) * t
00232                 + Q8) * t
00233                + Q7) * t
00234               + Q6) * t
00235               + Q5) * t
00236              + Q4) * t
00237             + Q3) * t
00238            + Q2) * t
00239           + Q1) * t
00240        + Q0;
00241       r = p / q;
00242       if (sign & 0x80000000)
00243        r = pimacosr4375 - r;
00244       else
00245        r = acosr4375 + r;
00246       return r;
00247     }
00248   else if (ix < 0x3ffe4000) /* |x| < 0.625 */
00249     {
00250       t = u.value - 0.5625L;
00251       p = ((((((((((rS10 * t
00252                   + rS9) * t
00253                  + rS8) * t
00254                 + rS7) * t
00255                + rS6) * t
00256               + rS5) * t
00257               + rS4) * t
00258              + rS3) * t
00259             + rS2) * t
00260            + rS1) * t
00261           + rS0) * t;
00262 
00263       q = (((((((((t
00264                  + sS9) * t
00265                 + sS8) * t
00266                + sS7) * t
00267               + sS6) * t
00268               + sS5) * t
00269              + sS4) * t
00270             + sS3) * t
00271            + sS2) * t
00272           + sS1) * t
00273        + sS0;
00274       if (sign & 0x80000000)
00275        r = pimacosr5625 - p / q;
00276       else
00277        r = acosr5625 + p / q;
00278       return r;
00279     }
00280   else
00281     {                       /* |x| >= .625 */
00282       z = (one - u.value) * 0.5;
00283       s = __ieee754_sqrtl (z);
00284       /* Compute an extended precision square root from
00285         the Newton iteration  s -> 0.5 * (s + z / s).
00286          The change w from s to the improved value is
00287            w = 0.5 * (s + z / s) - s  = (s^2 + z)/2s - s = (z - s^2)/2s.
00288           Express s = f1 + f2 where f1 * f1 is exactly representable.
00289          w = (z - s^2)/2s = (z - f1^2 - 2 f1 f2 - f2^2)/2s .
00290           s + w has extended precision.  */
00291       u.value = s;
00292       u.parts32.w2 = 0;
00293       u.parts32.w3 = 0;
00294       f2 = s - u.value;
00295       w = z - u.value * u.value;
00296       w = w - 2.0 * u.value * f2;
00297       w = w - f2 * f2;
00298       w = w / (2.0 * s);
00299       /* Arcsine of s.  */
00300       p = (((((((((pS9 * z
00301                  + pS8) * z
00302                 + pS7) * z
00303                + pS6) * z
00304               + pS5) * z
00305               + pS4) * z
00306              + pS3) * z
00307             + pS2) * z
00308            + pS1) * z
00309           + pS0) * z;
00310       q = (((((((( z
00311                  + qS8) * z
00312                + qS7) * z
00313               + qS6) * z
00314               + qS5) * z
00315              + qS4) * z
00316             + qS3) * z
00317            + qS2) * z
00318           + qS1) * z
00319        + qS0;
00320       r = s + (w + s * p / q);
00321 
00322       if (sign & 0x80000000)
00323        w = pio2_hi + (pio2_lo - r);
00324       else
00325        w = r;
00326       return 2.0 * w;
00327     }
00328 }