sysdeps/ieee754/ldbl-128ibm/e_acosl.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
  18    the 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 /* __ieee754_acosl(x)
  35  * Method :
  36  *      acos(x)  = pi/2 - asin(x)
  37  *      acos(-x) = pi/2 + asin(x)
  38  * For |x| <= 0.375
  39  *      acos(x) = pi/2 - asin(x)
  40  * Between .375 and .5 the approximation is
  41  *      acos(0.4375 + x) = acos(0.4375) + x P(x) / Q(x)
  42  * Between .5 and .625 the approximation is
  43  *      acos(0.5625 + x) = acos(0.5625) + x rS(x) / sS(x)
  44  * For x > 0.625,
  45  *      acos(x) = 2 asin(sqrt((1-x)/2))
  46  *      computed with an extended precision square root in the leading term.
  47  * For x < -0.625
  48  *      acos(x) = pi - 2 asin(sqrt((1-|x|)/2))
  49  *
  50  * Special cases:
  51  *      if x is NaN, return x itself;
  52  *      if |x|>1, return NaN with invalid signal.
  53  *
  54  * Functions needed: __ieee754_sqrtl.
  55  */
  56
  57 #include "math.h"
  58 #include "math_private.h"
  59
  60 #ifdef __STDC__
  61 static const long double
  62 #else
  63 static long double
  64 #endif
  65   one = 1.0L,
  66   pio2_hi = 1.5707963267948966192313216916397514420986L,
  67   pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
  68
  69   /* acos(0.5625 + x) = acos(0.5625) + x rS(x) / sS(x)
  70      -0.0625 <= x <= 0.0625
  71      peak relative error 3.3e-35  */
  72
  73   rS0 =  5.619049346208901520945464704848780243887E0L,
  74   rS1 = -4.460504162777731472539175700169871920352E1L,
  75   rS2 =  1.317669505315409261479577040530751477488E2L,
  76   rS3 = -1.626532582423661989632442410808596009227E2L,
  77   rS4 =  3.144806644195158614904369445440583873264E1L,
  78   rS5 =  9.806674443470740708765165604769099559553E1L,
  79   rS6 = -5.708468492052010816555762842394927806920E1L,
  80   rS7 = -1.396540499232262112248553357962639431922E1L,
  81   rS8 =  1.126243289311910363001762058295832610344E1L,
  82   rS9 =  4.956179821329901954211277873774472383512E-1L,
  83   rS10 = -3.313227657082367169241333738391762525780E-1L,
  84
  85   sS0 = -4.645814742084009935700221277307007679325E0L,
  86   sS1 =  3.879074822457694323970438316317961918430E1L,
  87   sS2 = -1.221986588013474694623973554726201001066E2L,
  88   sS3 =  1.658821150347718105012079876756201905822E2L,
  89   sS4 = -4.804379630977558197953176474426239748977E1L,
  90   sS5 = -1.004296417397316948114344573811562952793E2L,
  91   sS6 =  7.530281592861320234941101403870010111138E1L,
  92   sS7 =  1.270735595411673647119592092304357226607E1L,
  93   sS8 = -1.815144839646376500705105967064792930282E1L,
  94   sS9 = -7.821597334910963922204235247786840828217E-2L,
  95   /* 1.000000000000000000000000000000000000000E0 */
  96
  97   acosr5625 = 9.7338991014954640492751132535550279812151E-1L,
  98   pimacosr5625 = 2.1682027434402468335351320579240000860757E0L,
  99
 100   /* acos(0.4375 + x) = acos(0.4375) + x rS(x) / sS(x)
 101      -0.0625 <= x <= 0.0625
 102      peak relative error 2.1e-35  */
 103
 104   P0 =  2.177690192235413635229046633751390484892E0L,
 105   P1 = -2.848698225706605746657192566166142909573E1L,
 106   P2 =  1.040076477655245590871244795403659880304E2L,
 107   P3 = -1.400087608918906358323551402881238180553E2L,
 108   P4 =  2.221047917671449176051896400503615543757E1L,
 109   P5 =  9.643714856395587663736110523917499638702E1L,
 110   P6 = -5.158406639829833829027457284942389079196E1L,
 111   P7 = -1.578651828337585944715290382181219741813E1L,
 112   P8 =  1.093632715903802870546857764647931045906E1L,
 113   P9 =  5.448925479898460003048760932274085300103E-1L,
 114   P10 = -3.315886001095605268470690485170092986337E-1L,
 115   Q0 = -1.958219113487162405143608843774587557016E0L,
 116   Q1 =  2.614577866876185080678907676023269360520E1L,
 117   Q2 = -9.990858606464150981009763389881793660938E1L,
 118   Q3 =  1.443958741356995763628660823395334281596E2L,
 119   Q4 = -3.206441012484232867657763518369723873129E1L,
 120   Q5 = -1.048560885341833443564920145642588991492E2L,
 121   Q6 =  6.745883931909770880159915641984874746358E1L,
 122   Q7 =  1.806809656342804436118449982647641392951E1L,
 123   Q8 = -1.770150690652438294290020775359580915464E1L,
 124   Q9 = -5.659156469628629327045433069052560211164E-1L,
 125   /* 1.000000000000000000000000000000000000000E0 */
 126
 127   acosr4375 = 1.1179797320499710475919903296900511518755E0L,
 128   pimacosr4375 = 2.0236129215398221908706530535894517323217E0L,
 129
 130   /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
 131      0 <= x <= 0.5
 132      peak relative error 1.9e-35  */
 133   pS0 = -8.358099012470680544198472400254596543711E2L,
 134   pS1 =  3.674973957689619490312782828051860366493E3L,
 135   pS2 = -6.730729094812979665807581609853656623219E3L,
 136   pS3 =  6.643843795209060298375552684423454077633E3L,
 137   pS4 = -3.817341990928606692235481812252049415993E3L,
 138   pS5 =  1.284635388402653715636722822195716476156E3L,
 139   pS6 = -2.410736125231549204856567737329112037867E2L,
 140   pS7 =  2.219191969382402856557594215833622156220E1L,
 141   pS8 = -7.249056260830627156600112195061001036533E-1L,
 142   pS9 =  1.055923570937755300061509030361395604448E-3L,
 143
 144   qS0 = -5.014859407482408326519083440151745519205E3L,
 145   qS1 =  2.430653047950480068881028451580393430537E4L,
 146   qS2 = -4.997904737193653607449250593976069726962E4L,
 147   qS3 =  5.675712336110456923807959930107347511086E4L,
 148   qS4 = -3.881523118339661268482937768522572588022E4L,
 149   qS5 =  1.634202194895541569749717032234510811216E4L,
 150   qS6 = -4.151452662440709301601820849901296953752E3L,
 151   qS7 =  5.956050864057192019085175976175695342168E2L,
 152   qS8 = -4.175375777334867025769346564600396877176E1L;
 153   /* 1.000000000000000000000000000000000000000E0 */
 154
 155 #ifdef __STDC__
 156 long double
 157 __ieee754_acosl (long double x)
 158 #else
 159 long double
 160 __ieee754_acosl (x)
 161      long double x;
 162 #endif
 163 {
 164   long double z, r, w, p, q, s, t, f2;
 165   int32_t ix, sign;
 166   ieee854_long_double_shape_type u;
 167
 168   u.value = x;
 169   sign = u.parts32.w0;
 170   ix = sign & 0x7fffffff;
 171   u.parts32.w0 = ix;            /* |x| */
 172   if (ix >= 0x3ff00000)         /* |x| >= 1 */
 173     {
 174       if (ix == 0x3ff00000
 175           && (u.parts32.w1 | (u.parts32.w2&0x7fffffff) | u.parts32.w3) == 0)
 176         {                       /* |x| == 1 */
 177           if ((sign & 0x80000000) == 0)
 178             return 0.0;         /* acos(1) = 0  */
 179           else
 180             return (2.0 * pio2_hi) + (2.0 * pio2_lo);   /* acos(-1)= pi */
 181         }
 182       return (x - x) / (x - x); /* acos(|x| > 1) is NaN */
 183     }
 184   else if (ix < 0x3fe00000)     /* |x| < 0.5 */
 185     {
 186       if (ix < 0x3c600000)      /* |x| < 2**-57 */
 187         return pio2_hi + pio2_lo;
 188       if (ix < 0x3fde0000)      /* |x| < .4375 */
 189         {
 190           /* Arcsine of x.  */
 191           z = x * x;
 192           p = (((((((((pS9 * z
 193                        + pS8) * z
 194                       + pS7) * z
 195                      + pS6) * z
 196                     + pS5) * z
 197                    + pS4) * z
 198                   + pS3) * z
 199                  + pS2) * z
 200                 + pS1) * z
 201                + pS0) * z;
 202           q = (((((((( z
 203                        + qS8) * z
 204                      + qS7) * z
 205                     + qS6) * z
 206                    + qS5) * z
 207                   + qS4) * z
 208                  + qS3) * z
 209                 + qS2) * z
 210                + qS1) * z
 211             + qS0;
 212           r = x + x * p / q;
 213           z = pio2_hi - (r - pio2_lo);
 214           return z;
 215         }
 216       /* .4375 <= |x| < .5 */
 217       t = u.value - 0.4375L;
 218       p = ((((((((((P10 * t
 219                     + P9) * t
 220                    + P8) * t
 221                   + P7) * t
 222                  + P6) * t
 223                 + P5) * t
 224                + P4) * t
 225               + P3) * t
 226              + P2) * t
 227             + P1) * t
 228            + P0) * t;
 229
 230       q = (((((((((t
 231                    + Q9) * t
 232                   + Q8) * t
 233                  + Q7) * t
 234                 + Q6) * t
 235                + Q5) * t
 236               + Q4) * t
 237              + Q3) * t
 238             + Q2) * t
 239            + Q1) * t
 240         + Q0;
 241       r = p / q;
 242       if (sign & 0x80000000)
 243         r = pimacosr4375 - r;
 244       else
 245         r = acosr4375 + r;
 246       return r;
 247     }
 248   else if (ix < 0x3fe40000)     /* |x| < 0.625 */
 249     {
 250       t = u.value - 0.5625L;
 251       p = ((((((((((rS10 * t
 252                     + rS9) * t
 253                    + rS8) * t
 254                   + rS7) * t
 255                  + rS6) * t
 256                 + rS5) * t
 257                + rS4) * t
 258               + rS3) * t
 259              + rS2) * t
 260             + rS1) * t
 261            + rS0) * t;
 262
 263       q = (((((((((t
 264                    + sS9) * t
 265                   + sS8) * t
 266                  + sS7) * t
 267                 + sS6) * t
 268                + sS5) * t
 269               + sS4) * t
 270              + sS3) * t
 271             + sS2) * t
 272            + sS1) * t
 273         + sS0;
 274       if (sign & 0x80000000)
 275         r = pimacosr5625 - p / q;
 276       else
 277         r = acosr5625 + p / q;
 278       return r;
 279     }
 280   else
 281     {                           /* |x| >= .625 */
 282       z = (one - u.value) * 0.5;
 283       s = __ieee754_sqrtl (z);
 284       /* Compute an extended precision square root from
 285          the Newton iteration  s -> 0.5 * (s + z / s).
 286          The change w from s to the improved value is
 287             w = 0.5 * (s + z / s) - s  = (s^2 + z)/2s - s = (z - s^2)/2s.
 288           Express s = f1 + f2 where f1 * f1 is exactly representable.
 289           w = (z - s^2)/2s = (z - f1^2 - 2 f1 f2 - f2^2)/2s .
 290           s + w has extended precision.  */
 291       u.value = s;
 292       u.parts32.w2 = 0;
 293       u.parts32.w3 = 0;
 294       f2 = s - u.value;
 295       w = z - u.value * u.value;
 296       w = w - 2.0 * u.value * f2;
 297       w = w - f2 * f2;
 298       w = w / (2.0 * s);
 299       /* Arcsine of s.  */
 300       p = (((((((((pS9 * z
 301                    + pS8) * z
 302                   + pS7) * z
 303                  + pS6) * z
 304                 + pS5) * z
 305                + pS4) * z
 306               + pS3) * z
 307              + pS2) * z
 308             + pS1) * z
 309            + pS0) * z;
 310       q = (((((((( z
 311                    + qS8) * z
 312                  + qS7) * z
 313                 + qS6) * z
 314                + qS5) * z
 315               + qS4) * z
 316              + qS3) * z
 317             + qS2) * z
 318            + qS1) * z
 319         + qS0;
 320       r = s + (w + s * p / q);
 321
 322       if (sign & 0x80000000)
 323         w = pio2_hi + (pio2_lo - r);
 324       else
 325         w = r;
 326       return 2.0 * w;
 327     }
 328 }