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