sysdeps/ieee754/ldbl-96/e_j0l.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 /* Long double expansions are
  13   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
  14   and are incorporated herein by permission of the author.  The author
  15   reserves the right to distribute this material elsewhere under different
  16   copying permissions.  These modifications are distributed here under
  17   the following terms:
  18
  19     This library is free software; you can redistribute it and/or
  20     modify it under the terms of the GNU Lesser General Public
  21     License as published by the Free Software Foundation; either
  22     version 2.1 of the License, or (at your option) any later version.
  23
  24     This library is distributed in the hope that it will be useful,
  25     but WITHOUT ANY WARRANTY; without even the implied warranty of
  26     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  27     Lesser General Public License for more details.
  28
  29     You should have received a copy of the GNU Lesser General Public
  30     License along with this library; if not, write to the Free Software
  31     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
  32
  33 /* __ieee754_j0(x), __ieee754_y0(x)
  34  * Bessel function of the first and second kinds of order zero.
  35  * Method -- j0(x):
  36  *      1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
  37  *      2. Reduce x to |x| since j0(x)=j0(-x),  and
  38  *         for x in (0,2)
  39  *              j0(x) = 1 - z/4 + z^2*R0/S0,  where z = x*x;
  40  *         for x in (2,inf)
  41  *              j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
  42  *         where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
  43  *         as follow:
  44  *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
  45  *                      = 1/sqrt(2) * (cos(x) + sin(x))
  46  *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
  47  *                      = 1/sqrt(2) * (sin(x) - cos(x))
  48  *         (To avoid cancellation, use
  49  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  50  *          to compute the worse one.)
  51  *
  52  *      3 Special cases
  53  *              j0(nan)= nan
  54  *              j0(0) = 1
  55  *              j0(inf) = 0
  56  *
  57  * Method -- y0(x):
  58  *      1. For x<2.
  59  *         Since
  60  *              y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
  61  *         therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
  62  *         We use the following function to approximate y0,
  63  *              y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
  64  *
  65  *         Note: For tiny x, U/V = u0 and j0(x)~1, hence
  66  *              y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
  67  *      2. For x>=2.
  68  *              y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
  69  *         where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
  70  *         by the method mentioned above.
  71  *      3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
  72  */
  73
  74 #include "math.h"
  75 #include "math_private.h"
  76
  77 #ifdef __STDC__
  78 static long double pzero (long double), qzero (long double);
  79 #else
  80 static long double pzero (), qzero ();
  81 #endif
  82
  83 #ifdef __STDC__
  84 static const long double
  85 #else
  86 static long double
  87 #endif
  88   huge = 1e4930L,
  89   one = 1.0L,
  90   invsqrtpi = 5.6418958354775628694807945156077258584405e-1L,
  91   tpi = 6.3661977236758134307553505349005744813784e-1L,
  92
  93   /* J0(x) = 1 - x^2 / 4 + x^4 R0(x^2) / S0(x^2)
  94      0 <= x <= 2
  95      peak relative error 1.41e-22 */
  96   R[5] = {
  97   4.287176872744686992880841716723478740566E7L,
  98   -6.652058897474241627570911531740907185772E5L,
  99   7.011848381719789863458364584613651091175E3L,
 100   -3.168040850193372408702135490809516253693E1L,
 101   6.030778552661102450545394348845599300939E-2L,
 102 },
 103
 104  S[4] = {
 105    2.743793198556599677955266341699130654342E9L,
 106    3.364330079384816249840086842058954076201E7L,
 107    1.924119649412510777584684927494642526573E5L,
 108    6.239282256012734914211715620088714856494E2L,
 109    /*   1.000000000000000000000000000000000000000E0L,*/
 110 };
 111
 112 #ifdef __STDC__
 113 static const long double zero = 0.0;
 114 #else
 115 static long double zero = 0.0;
 116 #endif
 117
 118 #ifdef __STDC__
 119 long double
 120 __ieee754_j0l (long double x)
 121 #else
 122 long double
 123 __ieee754_j0l (x)
 124      long double x;
 125 #endif
 126 {
 127   long double z, s, c, ss, cc, r, u, v;
 128   int32_t ix;
 129   u_int32_t se, i0, i1;
 130
 131   GET_LDOUBLE_WORDS (se, i0, i1, x);
 132   ix = se & 0x7fff;
 133   if (ix >= 0x7fff)
 134     return one / (x * x);
 135   x = fabsl (x);
 136   if (ix >= 0x4000)             /* |x| >= 2.0 */
 137     {
 138       __sincosl (x, &s, &c);
 139       ss = s - c;
 140       cc = s + c;
 141       if (ix < 0x7ffe)
 142         {                       /* make sure x+x not overflow */
 143           z = -__cosl (x + x);
 144           if ((s * c) < zero)
 145             cc = z / ss;
 146           else
 147             ss = z / cc;
 148         }
 149       /*
 150        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 151        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 152        */
 153       if (ix > 0x4080)  /* 2^129 */
 154         z = (invsqrtpi * cc) / __ieee754_sqrtl (x);
 155       else
 156         {
 157           u = pzero (x);
 158           v = qzero (x);
 159           z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (x);
 160         }
 161       return z;
 162     }
 163   if (ix < 0x3fef) /* |x| < 2**-16 */
 164     {
 165       if (huge + x > one)
 166         {                       /* raise inexact if x != 0 */
 167           if (ix < 0x3fde) /* |x| < 2^-33 */
 168             return one;
 169           else
 170             return one - 0.25 * x * x;
 171         }
 172     }
 173   z = x * x;
 174   r = z * (R[0] + z * (R[1] + z * (R[2] + z * (R[3] + z * R[4]))));
 175   s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z)));
 176   if (ix < 0x3fff)
 177     {                           /* |x| < 1.00 */
 178       return (one - 0.25 * z + z * (r / s));
 179     }
 180   else
 181     {
 182       u = 0.5 * x;
 183       return ((one + u) * (one - u) + z * (r / s));
 184     }
 185 }
 186
 187
 188 /* y0(x) = 2/pi ln(x) J0(x) + U(x^2)/V(x^2)
 189    0 < x <= 2
 190    peak relative error 1.7e-21 */
 191 #ifdef __STDC__
 192 static const long double
 193 #else
 194 static long double
 195 #endif
 196 U[6] = {
 197   -1.054912306975785573710813351985351350861E10L,
 198   2.520192609749295139432773849576523636127E10L,
 199   -1.856426071075602001239955451329519093395E9L,
 200   4.079209129698891442683267466276785956784E7L,
 201   -3.440684087134286610316661166492641011539E5L,
 202   1.005524356159130626192144663414848383774E3L,
 203 };
 204 #ifdef __STDC__
 205 static const long double
 206 #else
 207 static long double
 208 #endif
 209 V[5] = {
 210   1.429337283720789610137291929228082613676E11L,
 211   2.492593075325119157558811370165695013002E9L,
 212   2.186077620785925464237324417623665138376E7L,
 213   1.238407896366385175196515057064384929222E5L,
 214   4.693924035211032457494368947123233101664E2L,
 215   /*  1.000000000000000000000000000000000000000E0L */
 216 };
 217
 218 #ifdef __STDC__
 219 long double
 220 __ieee754_y0l (long double x)
 221 #else
 222 long double
 223 __ieee754_y0l (x)
 224      long double x;
 225 #endif
 226 {
 227   long double z, s, c, ss, cc, u, v;
 228   int32_t ix;
 229   u_int32_t se, i0, i1;
 230
 231   GET_LDOUBLE_WORDS (se, i0, i1, x);
 232   ix = se & 0x7fff;
 233   /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
 234   if (se & 0x8000)
 235     return zero / (zero * x);
 236   if (ix >= 0x7fff)
 237     return one / (x + x * x);
 238   if ((i0 | i1) == 0)
 239     return -HUGE_VALL + x;  /* -inf and overflow exception.  */
 240   if (ix >= 0x4000)
 241     {                           /* |x| >= 2.0 */
 242
 243       /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
 244        * where x0 = x-pi/4
 245        *      Better formula:
 246        *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
 247        *                      =  1/sqrt(2) * (sin(x) + cos(x))
 248        *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 249        *                      =  1/sqrt(2) * (sin(x) - cos(x))
 250        * To avoid cancellation, use
 251        *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 252        * to compute the worse one.
 253        */
 254       __sincosl (x, &s, &c);
 255       ss = s - c;
 256       cc = s + c;
 257       /*
 258        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 259        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 260        */
 261       if (ix < 0x7ffe)
 262         {                       /* make sure x+x not overflow */
 263           z = -__cosl (x + x);
 264           if ((s * c) < zero)
 265             cc = z / ss;
 266           else
 267             ss = z / cc;
 268         }
 269       if (ix > 0x4080)  /* 1e39 */
 270         z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
 271       else
 272         {
 273           u = pzero (x);
 274           v = qzero (x);
 275           z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x);
 276         }
 277       return z;
 278     }
 279   if (ix <= 0x3fde) /* x < 2^-33 */
 280     {
 281       z = -7.380429510868722527629822444004602747322E-2L
 282         + tpi * __ieee754_logl (x);
 283       return z;
 284     }
 285   z = x * x;
 286   u = U[0] + z * (U[1] + z * (U[2] + z * (U[3] + z * (U[4] + z * U[5]))));
 287   v = V[0] + z * (V[1] + z * (V[2] + z * (V[3] + z * (V[4] + z))));
 288   return (u / v + tpi * (__ieee754_j0l (x) * __ieee754_logl (x)));
 289 }
 290
 291 /* The asymptotic expansions of pzero is
 292  *      1 - 9/128 s^2 + 11025/98304 s^4 - ...,  where s = 1/x.
 293  * For x >= 2, We approximate pzero by
 294  *      pzero(x) = 1 + s^2 R(s^2) / S(s^2)
 295  */
 296 #ifdef __STDC__
 297 static const long double pR8[7] = {
 298 #else
 299 static long double pR8[7] = {
 300 #endif
 301   /* 8 <= x <= inf
 302      Peak relative error 4.62 */
 303   -4.094398895124198016684337960227780260127E-9L,
 304   -8.929643669432412640061946338524096893089E-7L,
 305   -6.281267456906136703868258380673108109256E-5L,
 306   -1.736902783620362966354814353559382399665E-3L,
 307   -1.831506216290984960532230842266070146847E-2L,
 308   -5.827178869301452892963280214772398135283E-2L,
 309   -2.087563267939546435460286895807046616992E-2L,
 310 };
 311 #ifdef __STDC__
 312 static const long double pS8[6] = {
 313 #else
 314 static long double pS8[6] = {
 315 #endif
 316   5.823145095287749230197031108839653988393E-8L,
 317   1.279281986035060320477759999428992730280E-5L,
 318   9.132668954726626677174825517150228961304E-4L,
 319   2.606019379433060585351880541545146252534E-2L,
 320   2.956262215119520464228467583516287175244E-1L,
 321   1.149498145388256448535563278632697465675E0L,
 322   /* 1.000000000000000000000000000000000000000E0L, */
 323 };
 324
 325 #ifdef __STDC__
 326 static const long double pR5[7] = {
 327 #else
 328 static long double pR5[7] = {
 329 #endif
 330   /* 4.54541015625 <= x <= 8
 331      Peak relative error 6.51E-22 */
 332   -2.041226787870240954326915847282179737987E-7L,
 333   -2.255373879859413325570636768224534428156E-5L,
 334   -7.957485746440825353553537274569102059990E-4L,
 335   -1.093205102486816696940149222095559439425E-2L,
 336   -5.657957849316537477657603125260701114646E-2L,
 337   -8.641175552716402616180994954177818461588E-2L,
 338   -1.354654710097134007437166939230619726157E-2L,
 339 };
 340 #ifdef __STDC__
 341 static const long double pS5[6] = {
 342 #else
 343 static long double pS5[6] = {
 344 #endif
 345   2.903078099681108697057258628212823545290E-6L,
 346   3.253948449946735405975737677123673867321E-4L,
 347   1.181269751723085006534147920481582279979E-2L,
 348   1.719212057790143888884745200257619469363E-1L,
 349   1.006306498779212467670654535430694221924E0L,
 350   2.069568808688074324555596301126375951502E0L,
 351   /* 1.000000000000000000000000000000000000000E0L, */
 352 };
 353
 354 #ifdef __STDC__
 355 static const long double pR3[7] = {
 356 #else
 357 static long double pR3[7] = {
 358 #endif
 359   /* 2.85711669921875 <= x <= 4.54541015625
 360      peak relative error 5.25e-21 */
 361   -5.755732156848468345557663552240816066802E-6L,
 362   -3.703675625855715998827966962258113034767E-4L,
 363   -7.390893350679637611641350096842846433236E-3L,
 364   -5.571922144490038765024591058478043873253E-2L,
 365   -1.531290690378157869291151002472627396088E-1L,
 366   -1.193350853469302941921647487062620011042E-1L,
 367   -8.567802507331578894302991505331963782905E-3L,
 368 };
 369 #ifdef __STDC__
 370 static const long double pS3[6] = {
 371 #else
 372 static long double pS3[6] = {
 373 #endif
 374   8.185931139070086158103309281525036712419E-5L,
 375   5.398016943778891093520574483111255476787E-3L,
 376   1.130589193590489566669164765853409621081E-1L,
 377   9.358652328786413274673192987670237145071E-1L,
 378   3.091711512598349056276917907005098085273E0L,
 379   3.594602474737921977972586821673124231111E0L,
 380   /* 1.000000000000000000000000000000000000000E0L, */
 381 };
 382
 383 #ifdef __STDC__
 384 static const long double pR2[7] = {
 385 #else
 386 static long double pR2[7] = {
 387 #endif
 388   /* 2 <= x <= 2.85711669921875
 389      peak relative error 2.64e-21 */
 390   -1.219525235804532014243621104365384992623E-4L,
 391   -4.838597135805578919601088680065298763049E-3L,
 392   -5.732223181683569266223306197751407418301E-2L,
 393   -2.472947430526425064982909699406646503758E-1L,
 394   -3.753373645974077960207588073975976327695E-1L,
 395   -1.556241316844728872406672349347137975495E-1L,
 396   -5.355423239526452209595316733635519506958E-3L,
 397 };
 398 #ifdef __STDC__
 399 static const long double pS2[6] = {
 400 #else
 401 static long double pS2[6] = {
 402 #endif
 403   1.734442793664291412489066256138894953823E-3L,
 404   7.158111826468626405416300895617986926008E-2L,
 405   9.153839713992138340197264669867993552641E-1L,
 406   4.539209519433011393525841956702487797582E0L,
 407   8.868932430625331650266067101752626253644E0L,
 408   6.067161890196324146320763844772857713502E0L,
 409   /* 1.000000000000000000000000000000000000000E0L, */
 410 };
 411
 412 #ifdef __STDC__
 413 static long double
 414 pzero (long double x)
 415 #else
 416 static long double
 417 pzero (x)
 418      long double x;
 419 #endif
 420 {
 421 #ifdef __STDC__
 422   const long double *p, *q;
 423 #else
 424   long double *p, *q;
 425 #endif
 426   long double z, r, s;
 427   int32_t ix;
 428   u_int32_t se, i0, i1;
 429
 430   GET_LDOUBLE_WORDS (se, i0, i1, x);
 431   ix = se & 0x7fff;
 432   if (ix >= 0x4002)
 433     {
 434       p = pR8;
 435       q = pS8;
 436     }                           /* x >= 8 */
 437   else
 438     {
 439       i1 = (ix << 16) | (i0 >> 16);
 440       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 441         {
 442           p = pR5;
 443           q = pS5;
 444         }
 445       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 446         {
 447           p = pR3;
 448           q = pS3;
 449         }
 450       else if (ix >= 0x4000)    /* x better be >= 2 */
 451         {
 452           p = pR2;
 453           q = pS2;
 454         }
 455     }
 456   z = one / (x * x);
 457   r =
 458     p[0] + z * (p[1] +
 459                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 460   s =
 461     q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
 462   return (one + z * r / s);
 463 }
 464
 465
 466 /* For x >= 8, the asymptotic expansions of qzero is
 467  *      -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
 468  * We approximate qzero by
 469  *      qzero(x) = s*(-.125 + R(s^2) / S(s^2))
 470  */
 471 #ifdef __STDC__
 472 static const long double qR8[7] = {
 473 #else
 474 static long double qR8[7] = {
 475 #endif
 476   /* 8 <= x <= inf
 477      peak relative error 2.23e-21 */
 478   3.001267180483191397885272640777189348008E-10L,
 479   8.693186311430836495238494289942413810121E-8L,
 480   8.496875536711266039522937037850596580686E-6L,
 481   3.482702869915288984296602449543513958409E-4L,
 482   6.036378380706107692863811938221290851352E-3L,
 483   3.881970028476167836382607922840452192636E-2L,
 484   6.132191514516237371140841765561219149638E-2L,
 485 };
 486 #ifdef __STDC__
 487 static const long double qS8[7] = {
 488 #else
 489 static long double qS8[7] = {
 490 #endif
 491   4.097730123753051126914971174076227600212E-9L,
 492   1.199615869122646109596153392152131139306E-6L,
 493   1.196337580514532207793107149088168946451E-4L,
 494   5.099074440112045094341500497767181211104E-3L,
 495   9.577420799632372483249761659674764460583E-2L,
 496   7.385243015344292267061953461563695918646E-1L,
 497   1.917266424391428937962682301561699055943E0L,
 498   /* 1.000000000000000000000000000000000000000E0L, */
 499 };
 500
 501 #ifdef __STDC__
 502 static const long double qR5[7] = {
 503 #else
 504 static long double qR5[7] = {
 505 #endif
 506   /* 4.54541015625 <= x <= 8
 507      peak relative error 1.03e-21 */
 508   3.406256556438974327309660241748106352137E-8L,
 509   4.855492710552705436943630087976121021980E-6L,
 510   2.301011739663737780613356017352912281980E-4L,
 511   4.500470249273129953870234803596619899226E-3L,
 512   3.651376459725695502726921248173637054828E-2L,
 513   1.071578819056574524416060138514508609805E-1L,
 514   7.458950172851611673015774675225656063757E-2L,
 515 };
 516 #ifdef __STDC__
 517 static const long double qS5[7] = {
 518 #else
 519 static long double qS5[7] = {
 520 #endif
 521   4.650675622764245276538207123618745150785E-7L,
 522   6.773573292521412265840260065635377164455E-5L,
 523   3.340711249876192721980146877577806687714E-3L,
 524   7.036218046856839214741678375536970613501E-2L,
 525   6.569599559163872573895171876511377891143E-1L,
 526   2.557525022583599204591036677199171155186E0L,
 527   3.457237396120935674982927714210361269133E0L,
 528   /* 1.000000000000000000000000000000000000000E0L,*/
 529 };
 530
 531 #ifdef __STDC__
 532 static const long double qR3[7] = {
 533 #else
 534 static long double qR3[7] = {
 535 #endif
 536   /* 2.85711669921875 <= x <= 4.54541015625
 537      peak relative error 5.24e-21 */
 538   1.749459596550816915639829017724249805242E-6L,
 539   1.446252487543383683621692672078376929437E-4L,
 540   3.842084087362410664036704812125005761859E-3L,
 541   4.066369994699462547896426554180954233581E-2L,
 542   1.721093619117980251295234795188992722447E-1L,
 543   2.538595333972857367655146949093055405072E-1L,
 544   8.560591367256769038905328596020118877936E-2L,
 545 };
 546 #ifdef __STDC__
 547 static const long double qS3[7] = {
 548 #else
 549 static long double qS3[7] = {
 550 #endif
 551   2.388596091707517488372313710647510488042E-5L,
 552   2.048679968058758616370095132104333998147E-3L,
 553   5.824663198201417760864458765259945181513E-2L,
 554   6.953906394693328750931617748038994763958E-1L,
 555   3.638186936390881159685868764832961092476E0L,
 556   7.900169524705757837298990558459547842607E0L,
 557   5.992718532451026507552820701127504582907E0L,
 558   /* 1.000000000000000000000000000000000000000E0L, */
 559 };
 560
 561 #ifdef __STDC__
 562 static const long double qR2[7] = {
 563 #else
 564 static long double qR2[7] = {
 565 #endif
 566   /* 2 <= x <= 2.85711669921875
 567      peak relative error 1.58e-21  */
 568   6.306524405520048545426928892276696949540E-5L,
 569   3.209606155709930950935893996591576624054E-3L,
 570   5.027828775702022732912321378866797059604E-2L,
 571   3.012705561838718956481911477587757845163E-1L,
 572   6.960544893905752937420734884995688523815E-1L,
 573   5.431871999743531634887107835372232030655E-1L,
 574   9.447736151202905471899259026430157211949E-2L,
 575 };
 576 #ifdef __STDC__
 577 static const long double qS2[7] = {
 578 #else
 579 static long double qS2[7] = {
 580 #endif
 581   8.610579901936193494609755345106129102676E-4L,
 582   4.649054352710496997203474853066665869047E-2L,
 583   8.104282924459837407218042945106320388339E-1L,
 584   5.807730930825886427048038146088828206852E0L,
 585   1.795310145936848873627710102199881642939E1L,
 586   2.281313316875375733663657188888110605044E1L,
 587   1.011242067883822301487154844458322200143E1L,
 588   /* 1.000000000000000000000000000000000000000E0L, */
 589 };
 590
 591 #ifdef __STDC__
 592 static long double
 593 qzero (long double x)
 594 #else
 595 static long double
 596 qzero (x)
 597      long double x;
 598 #endif
 599 {
 600 #ifdef __STDC__
 601   const long double *p, *q;
 602 #else
 603   long double *p, *q;
 604 #endif
 605   long double s, r, z;
 606   int32_t ix;
 607   u_int32_t se, i0, i1;
 608
 609   GET_LDOUBLE_WORDS (se, i0, i1, x);
 610   ix = se & 0x7fff;
 611   if (ix >= 0x4002)             /* x >= 8 */
 612     {
 613       p = qR8;
 614       q = qS8;
 615     }
 616   else
 617     {
 618       i1 = (ix << 16) | (i0 >> 16);
 619       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 620         {
 621           p = qR5;
 622           q = qS5;
 623         }
 624       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 625         {
 626           p = qR3;
 627           q = qS3;
 628         }
 629       else if (ix >= 0x4000)    /* x better be >= 2 */
 630         {
 631           p = qR2;
 632           q = qS2;
 633         }
 634     }
 635   z = one / (x * x);
 636   r =
 637     p[0] + z * (p[1] +
 638                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 639   s =
 640     q[0] + z * (q[1] +
 641                 z * (q[2] +
 642                      z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
 643   return (-.125 + z * r / s) / x;
 644 }