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   a = __builtin_fabsl (x);
 157   if (a == 1.0L)
 158     {
 159       if (x > 0.0L)
 160         return 0.0;             /* acos(1) = 0  */
 161       else
 162         return (2.0 * pio2_hi) + (2.0 * pio2_lo);       /* acos(-1)= pi */
 163     }
 164   else if (a > 1.0L)
 165     {
 166       return (x - x) / (x - x); /* acos(|x| > 1) is NaN */
 167     }
 168   if (a < 0.5L)
 169     {
 170       if (a < 6.938893903907228e-18L)   /* |x| < 2**-57 */
 171         return pio2_hi + pio2_lo;
 172       if (a < 0.4375L)
 173         {
 174           /* Arcsine of x.  */
 175           z = x * x;
 176           p = (((((((((pS9 * z
 177                        + pS8) * z
 178                       + pS7) * z
 179                      + pS6) * z
 180                     + pS5) * z
 181                    + pS4) * z
 182                   + pS3) * z
 183                  + pS2) * z
 184                 + pS1) * z
 185                + pS0) * z;
 186           q = (((((((( z
 187                        + qS8) * z
 188                      + qS7) * z
 189                     + qS6) * z
 190                    + qS5) * z
 191                   + qS4) * z
 192                  + qS3) * z
 193                 + qS2) * z
 194                + qS1) * z
 195             + qS0;
 196           r = x + x * p / q;
 197           z = pio2_hi - (r - pio2_lo);
 198           return z;
 199         }
 200       /* .4375 <= |x| < .5 */
 201       t = a - 0.4375L;
 202       p = ((((((((((P10 * t
 203                     + P9) * t
 204                    + P8) * t
 205                   + P7) * t
 206                  + P6) * t
 207                 + P5) * t
 208                + P4) * t
 209               + P3) * t
 210              + P2) * t
 211             + P1) * t
 212            + P0) * t;
 213
 214       q = (((((((((t
 215                    + Q9) * t
 216                   + Q8) * t
 217                  + Q7) * t
 218                 + Q6) * t
 219                + Q5) * t
 220               + Q4) * t
 221              + Q3) * t
 222             + Q2) * t
 223            + Q1) * t
 224         + Q0;
 225       r = p / q;
 226       if (x < 0.0L)
 227         r = pimacosr4375 - r;
 228       else
 229         r = acosr4375 + r;
 230       return r;
 231     }
 232   else if (a < 0.625L)
 233     {
 234       t = a - 0.5625L;
 235       p = ((((((((((rS10 * t
 236                     + rS9) * t
 237                    + rS8) * t
 238                   + rS7) * t
 239                  + rS6) * t
 240                 + rS5) * t
 241                + rS4) * t
 242               + rS3) * t
 243              + rS2) * t
 244             + rS1) * t
 245            + rS0) * t;
 246
 247       q = (((((((((t
 248                    + sS9) * t
 249                   + sS8) * t
 250                  + sS7) * t
 251                 + sS6) * t
 252                + sS5) * t
 253               + sS4) * t
 254              + sS3) * t
 255             + sS2) * t
 256            + sS1) * t
 257         + sS0;
 258       if (x < 0.0L)
 259         r = pimacosr5625 - p / q;
 260       else
 261         r = acosr5625 + p / q;
 262       return r;
 263     }
 264   else
 265     {                           /* |x| >= .625 */
 266       double shi, slo;
 267
 268       z = (one - a) * 0.5;
 269       s = __ieee754_sqrtl (z);
 270       /* Compute an extended precision square root from
 271          the Newton iteration  s -> 0.5 * (s + z / s).
 272          The change w from s to the improved value is
 273             w = 0.5 * (s + z / s) - s  = (s^2 + z)/2s - s = (z - s^2)/2s.
 274           Express s = f1 + f2 where f1 * f1 is exactly representable.
 275           w = (z - s^2)/2s = (z - f1^2 - 2 f1 f2 - f2^2)/2s .
 276           s + w has extended precision.  */
 277       ldbl_unpack (s, &shi, &slo);
 278       a = shi;
 279       f2 = slo;
 280       w = z - a * a;
 281       w = w - 2.0 * a * f2;
 282       w = w - f2 * f2;
 283       w = w / (2.0 * s);
 284       /* Arcsine of s.  */
 285       p = (((((((((pS9 * z
 286                    + pS8) * z
 287                   + pS7) * z
 288                  + pS6) * z
 289                 + pS5) * z
 290                + pS4) * z
 291               + pS3) * z
 292              + pS2) * z
 293             + pS1) * z
 294            + pS0) * z;
 295       q = (((((((( z
 296                    + qS8) * z
 297                  + qS7) * z
 298                 + qS6) * z
 299                + qS5) * z
 300               + qS4) * z
 301              + qS3) * z
 302             + qS2) * z
 303            + qS1) * z
 304         + qS0;
 305       r = s + (w + s * p / q);
 306
 307       if (x < 0.0L)
 308         w = pio2_hi + (pio2_lo - r);
 309       else
 310         w = r;
 311       return 2.0 * w;
 312     }
 313 }
 314 strong_alias (__ieee754_acosl, __acosl_finite)