sysdeps/ieee754/ldbl-96/e_j1l.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_j1(x), __ieee754_y1(x)
  34  * Bessel function of the first and second kinds of order zero.
  35  * Method -- j1(x):
  36  *      1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
  37  *      2. Reduce x to |x| since j1(x)=-j1(-x),  and
  38  *         for x in (0,2)
  39  *              j1(x) = x/2 + x*z*R0/S0,  where z = x*x;
  40  *         for x in (2,inf)
  41  *              j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
  42  *              y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
  43  *         where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
  44  *         as follow:
  45  *              cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
  46  *                      =  1/sqrt(2) * (sin(x) - cos(x))
  47  *              sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
  48  *                      = -1/sqrt(2) * (sin(x) + cos(x))
  49  *         (To avoid cancellation, use
  50  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  51  *          to compute the worse one.)
  52  *
  53  *      3 Special cases
  54  *              j1(nan)= nan
  55  *              j1(0) = 0
  56  *              j1(inf) = 0
  57  *
  58  * Method -- y1(x):
  59  *      1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
  60  *      2. For x<2.
  61  *         Since
  62  *              y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
  63  *         therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
  64  *         We use the following function to approximate y1,
  65  *              y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
  66  *         Note: For tiny x, 1/x dominate y1 and hence
  67  *              y1(tiny) = -2/pi/tiny
  68  *      3. For x>=2.
  69  *              y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
  70  *         where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
  71  *         by method mentioned above.
  72  */
  73
  74 #include "math.h"
  75 #include "math_private.h"
  76
  77 #ifdef __STDC__
  78 static long double pone (long double), qone (long double);
  79 #else
  80 static long double pone (), qone ();
  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   /* J1(x) = .5 x + x x^2 R(x^2) / S(x^2)
  94      0 <= x <= 2
  95      Peak relative error 4.5e-21 */
  96 R[5] = {
  97     -9.647406112428107954753770469290757756814E7L,
  98     2.686288565865230690166454005558203955564E6L,
  99     -3.689682683905671185891885948692283776081E4L,
 100     2.195031194229176602851429567792676658146E2L,
 101     -5.124499848728030297902028238597308971319E-1L,
 102 },
 103
 104   S[4] =
 105 {
 106   1.543584977988497274437410333029029035089E9L,
 107   2.133542369567701244002565983150952549520E7L,
 108   1.394077011298227346483732156167414670520E5L,
 109   5.252401789085732428842871556112108446506E2L,
 110   /*  1.000000000000000000000000000000000000000E0L, */
 111 };
 112
 113 #ifdef __STDC__
 114 static const long double zero = 0.0;
 115 #else
 116 static long double zero = 0.0;
 117 #endif
 118
 119
 120 #ifdef __STDC__
 121 long double
 122 __ieee754_j1l (long double x)
 123 #else
 124 long double
 125 __ieee754_j1l (x)
 126      long double x;
 127 #endif
 128 {
 129   long double z, c, r, s, ss, cc, u, v, y;
 130   int32_t ix;
 131   u_int32_t se, i0, i1;
 132
 133   GET_LDOUBLE_WORDS (se, i0, i1, x);
 134   ix = se & 0x7fff;
 135   if (ix >= 0x7fff)
 136     return one / x;
 137   y = fabsl (x);
 138   if (ix >= 0x4000)
 139     {                           /* |x| >= 2.0 */
 140       __sincosl (y, &s, &c);
 141       ss = -s - c;
 142       cc = s - c;
 143       if (ix < 0x7ffe)
 144         {                       /* make sure y+y not overflow */
 145           z = __cosl (y + y);
 146           if ((s * c) > zero)
 147             cc = z / ss;
 148           else
 149             ss = z / cc;
 150         }
 151       /*
 152        * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
 153        * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
 154        */
 155       if (ix > 0x4080)
 156         z = (invsqrtpi * cc) / __ieee754_sqrtl (y);
 157       else
 158         {
 159           u = pone (y);
 160           v = qone (y);
 161           z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (y);
 162         }
 163       if (se & 0x8000)
 164         return -z;
 165       else
 166         return z;
 167     }
 168   if (ix < 0x3fde) /* |x| < 2^-33 */
 169     {
 170       if (huge + x > one)
 171         return 0.5 * x;         /* inexact if x!=0 necessary */
 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   r *= x;
 177   return (x * 0.5 + r / s);
 178 }
 179
 180
 181 /* Y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + x R(x^2)
 182    0 <= x <= 2
 183    Peak relative error 2.3e-23 */
 184 #ifdef __STDC__
 185 static const long double U0[6] = {
 186 #else
 187 static long double U0[6] = {
 188 #endif
 189   -5.908077186259914699178903164682444848615E10L,
 190   1.546219327181478013495975514375773435962E10L,
 191   -6.438303331169223128870035584107053228235E8L,
 192   9.708540045657182600665968063824819371216E6L,
 193   -6.138043997084355564619377183564196265471E4L,
 194   1.418503228220927321096904291501161800215E2L,
 195 };
 196 #ifdef __STDC__
 197 static const long double V0[5] = {
 198 #else
 199 static long double V0[5] = {
 200 #endif
 201   3.013447341682896694781964795373783679861E11L,
 202   4.669546565705981649470005402243136124523E9L,
 203   3.595056091631351184676890179233695857260E7L,
 204   1.761554028569108722903944659933744317994E5L,
 205   5.668480419646516568875555062047234534863E2L,
 206   /*  1.000000000000000000000000000000000000000E0L, */
 207 };
 208
 209
 210 #ifdef __STDC__
 211 long double
 212 __ieee754_y1l (long double x)
 213 #else
 214 long double
 215 __ieee754_y1l (x)
 216      long double x;
 217 #endif
 218 {
 219   long double z, s, c, ss, cc, u, v;
 220   int32_t ix;
 221   u_int32_t se, i0, i1;
 222
 223   GET_LDOUBLE_WORDS (se, i0, i1, x);
 224   ix = se & 0x7fff;
 225   /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
 226   if (se & 0x8000)
 227     return zero / (zero * x);
 228   if (ix >= 0x7fff)
 229     return one / (x + x * x);
 230   if ((i0 | i1) == 0)
 231     return -HUGE_VALL + x;  /* -inf and overflow exception.  */
 232   if (ix >= 0x4000)
 233     {                           /* |x| >= 2.0 */
 234       __sincosl (x, &s, &c);
 235       ss = -s - c;
 236       cc = s - c;
 237       if (ix < 0x7fe00000)
 238         {                       /* make sure x+x not overflow */
 239           z = __cosl (x + x);
 240           if ((s * c) > zero)
 241             cc = z / ss;
 242           else
 243             ss = z / cc;
 244         }
 245       /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
 246        * where x0 = x-3pi/4
 247        *      Better formula:
 248        *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
 249        *                      =  1/sqrt(2) * (sin(x) - cos(x))
 250        *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 251        *                      = -1/sqrt(2) * (cos(x) + sin(x))
 252        * To avoid cancellation, use
 253        *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 254        * to compute the worse one.
 255        */
 256       if (ix > 0x4080)
 257         z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
 258       else
 259         {
 260           u = pone (x);
 261           v = qone (x);
 262           z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x);
 263         }
 264       return z;
 265     }
 266   if (ix <= 0x3fbe)
 267     {                           /* x < 2**-65 */
 268       return (-tpi / x);
 269     }
 270   z = x * x;
 271  u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * (U0[4] + z * U0[5]))));
 272  v = V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * (V0[4] + z))));
 273   return (x * (u / v) +
 274           tpi * (__ieee754_j1l (x) * __ieee754_logl (x) - one / x));
 275 }
 276
 277
 278 /* For x >= 8, the asymptotic expansions of pone is
 279  *      1 + 15/128 s^2 - 4725/2^15 s^4 - ...,   where s = 1/x.
 280  * We approximate pone by
 281  *      pone(x) = 1 + (R/S)
 282  */
 283
 284 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 285    P1(x) = 1 + z^2 R(z^2), z=1/x
 286    8 <= x <= inf  (0 <= z <= 0.125)
 287    Peak relative error 5.2e-22  */
 288
 289 #ifdef __STDC__
 290 static const long double pr8[7] = {
 291 #else
 292 static long double pr8[7] = {
 293 #endif
 294   8.402048819032978959298664869941375143163E-9L,
 295   1.813743245316438056192649247507255996036E-6L,
 296   1.260704554112906152344932388588243836276E-4L,
 297   3.439294839869103014614229832700986965110E-3L,
 298   3.576910849712074184504430254290179501209E-2L,
 299   1.131111483254318243139953003461511308672E-1L,
 300   4.480715825681029711521286449131671880953E-2L,
 301 };
 302 #ifdef __STDC__
 303 static const long double ps8[6] = {
 304 #else
 305 static long double ps8[6] = {
 306 #endif
 307   7.169748325574809484893888315707824924354E-8L,
 308   1.556549720596672576431813934184403614817E-5L,
 309   1.094540125521337139209062035774174565882E-3L,
 310   3.060978962596642798560894375281428805840E-2L,
 311   3.374146536087205506032643098619414507024E-1L,
 312   1.253830208588979001991901126393231302559E0L,
 313   /* 1.000000000000000000000000000000000000000E0L, */
 314 };
 315
 316 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 317    P1(x) = 1 + z^2 R(z^2), z=1/x
 318    4.54541015625 <= x <= 8
 319    Peak relative error 7.7e-22  */
 320 #ifdef __STDC__
 321 static const long double pr5[7] = {
 322 #else
 323 static long double pr5[7] = {
 324 #endif
 325   4.318486887948814529950980396300969247900E-7L,
 326   4.715341880798817230333360497524173929315E-5L,
 327   1.642719430496086618401091544113220340094E-3L,
 328   2.228688005300803935928733750456396149104E-2L,
 329   1.142773760804150921573259605730018327162E-1L,
 330   1.755576530055079253910829652698703791957E-1L,
 331   3.218803858282095929559165965353784980613E-2L,
 332 };
 333 #ifdef __STDC__
 334 static const long double ps5[6] = {
 335 #else
 336 static long double ps5[6] = {
 337 #endif
 338   3.685108812227721334719884358034713967557E-6L,
 339   4.069102509511177498808856515005792027639E-4L,
 340   1.449728676496155025507893322405597039816E-2L,
 341   2.058869213229520086582695850441194363103E-1L,
 342   1.164890985918737148968424972072751066553E0L,
 343   2.274776933457009446573027260373361586841E0L,
 344   /*  1.000000000000000000000000000000000000000E0L,*/
 345 };
 346
 347 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 348    P1(x) = 1 + z^2 R(z^2), z=1/x
 349    2.85711669921875 <= x <= 4.54541015625
 350    Peak relative error 6.5e-21  */
 351 #ifdef __STDC__
 352 static const long double pr3[7] = {
 353 #else
 354 static long double pr3[7] = {
 355 #endif
 356   1.265251153957366716825382654273326407972E-5L,
 357   8.031057269201324914127680782288352574567E-4L,
 358   1.581648121115028333661412169396282881035E-2L,
 359   1.179534658087796321928362981518645033967E-1L,
 360   3.227936912780465219246440724502790727866E-1L,
 361   2.559223765418386621748404398017602935764E-1L,
 362   2.277136933287817911091370397134882441046E-2L,
 363 };
 364 #ifdef __STDC__
 365 static const long double ps3[6] = {
 366 #else
 367 static long double ps3[6] = {
 368 #endif
 369   1.079681071833391818661952793568345057548E-4L,
 370   6.986017817100477138417481463810841529026E-3L,
 371   1.429403701146942509913198539100230540503E-1L,
 372   1.148392024337075609460312658938700765074E0L,
 373   3.643663015091248720208251490291968840882E0L,
 374   3.990702269032018282145100741746633960737E0L,
 375   /* 1.000000000000000000000000000000000000000E0L, */
 376 };
 377
 378 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 379    P1(x) = 1 + z^2 R(z^2), z=1/x
 380    2 <= x <= 2.85711669921875
 381    Peak relative error 3.5e-21  */
 382 #ifdef __STDC__
 383 static const long double pr2[7] = {
 384 #else
 385 static long double pr2[7] = {
 386 #endif
 387   2.795623248568412225239401141338714516445E-4L,
 388   1.092578168441856711925254839815430061135E-2L,
 389   1.278024620468953761154963591853679640560E-1L,
 390   5.469680473691500673112904286228351988583E-1L,
 391   8.313769490922351300461498619045639016059E-1L,
 392   3.544176317308370086415403567097130611468E-1L,
 393   1.604142674802373041247957048801599740644E-2L,
 394 };
 395 #ifdef __STDC__
 396 static const long double ps2[6] = {
 397 #else
 398 static long double ps2[6] = {
 399 #endif
 400   2.385605161555183386205027000675875235980E-3L,
 401   9.616778294482695283928617708206967248579E-2L,
 402   1.195215570959693572089824415393951258510E0L,
 403   5.718412857897054829999458736064922974662E0L,
 404   1.065626298505499086386584642761602177568E1L,
 405   6.809140730053382188468983548092322151791E0L,
 406  /* 1.000000000000000000000000000000000000000E0L, */
 407 };
 408
 409
 410 #ifdef __STDC__
 411 static long double
 412 pone (long double x)
 413 #else
 414 static long double
 415 pone (x)
 416      long double x;
 417 #endif
 418 {
 419 #ifdef __STDC__
 420   const long double *p, *q;
 421 #else
 422   long double *p, *q;
 423 #endif
 424   long double z, r, s;
 425   int32_t ix;
 426   u_int32_t se, i0, i1;
 427
 428   GET_LDOUBLE_WORDS (se, i0, i1, x);
 429   ix = se & 0x7fff;
 430   if (ix >= 0x4002) /* x >= 8 */
 431     {
 432       p = pr8;
 433       q = ps8;
 434     }
 435   else
 436     {
 437       i1 = (ix << 16) | (i0 >> 16);
 438       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 439         {
 440           p = pr5;
 441           q = ps5;
 442         }
 443       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 444         {
 445           p = pr3;
 446           q = ps3;
 447         }
 448       else if (ix >= 0x4000)    /* x better be >= 2 */
 449         {
 450           p = pr2;
 451           q = ps2;
 452         }
 453     }
 454   z = one / (x * x);
 455  r = p[0] + z * (p[1] +
 456                  z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 457  s = q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
 458   return one + z * r / s;
 459 }
 460
 461
 462 /* For x >= 8, the asymptotic expansions of qone is
 463  *      3/8 s - 105/1024 s^3 - ..., where s = 1/x.
 464  * We approximate pone by
 465  *      qone(x) = s*(0.375 + (R/S))
 466  */
 467
 468 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 469    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 470    8 <= x <= inf
 471    Peak relative error 8.3e-22 */
 472
 473 #ifdef __STDC__
 474 static const long double qr8[7] = {
 475 #else
 476 static long double qr8[7] = {
 477 #endif
 478   -5.691925079044209246015366919809404457380E-10L,
 479   -1.632587664706999307871963065396218379137E-7L,
 480   -1.577424682764651970003637263552027114600E-5L,
 481   -6.377627959241053914770158336842725291713E-4L,
 482   -1.087408516779972735197277149494929568768E-2L,
 483   -6.854943629378084419631926076882330494217E-2L,
 484   -1.055448290469180032312893377152490183203E-1L,
 485 };
 486 #ifdef __STDC__
 487 static const long double qs8[7] = {
 488 #else
 489 static long double qs8[7] = {
 490 #endif
 491   5.550982172325019811119223916998393907513E-9L,
 492   1.607188366646736068460131091130644192244E-6L,
 493   1.580792530091386496626494138334505893599E-4L,
 494   6.617859900815747303032860443855006056595E-3L,
 495   1.212840547336984859952597488863037659161E-1L,
 496   9.017885953937234900458186716154005541075E-1L,
 497   2.201114489712243262000939120146436167178E0L,
 498   /* 1.000000000000000000000000000000000000000E0L, */
 499 };
 500
 501 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 502    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 503    4.54541015625 <= x <= 8
 504    Peak relative error 4.1e-22 */
 505 #ifdef __STDC__
 506 static const long double qr5[7] = {
 507 #else
 508 static long double qr5[7] = {
 509 #endif
 510   -6.719134139179190546324213696633564965983E-8L,
 511   -9.467871458774950479909851595678622044140E-6L,
 512   -4.429341875348286176950914275723051452838E-4L,
 513   -8.539898021757342531563866270278505014487E-3L,
 514   -6.818691805848737010422337101409276287170E-2L,
 515   -1.964432669771684034858848142418228214855E-1L,
 516   -1.333896496989238600119596538299938520726E-1L,
 517 };
 518 #ifdef __STDC__
 519 static const long double qs5[7] = {
 520 #else
 521 static long double qs5[7] = {
 522 #endif
 523   6.552755584474634766937589285426911075101E-7L,
 524   9.410814032118155978663509073200494000589E-5L,
 525   4.561677087286518359461609153655021253238E-3L,
 526   9.397742096177905170800336715661091535805E-2L,
 527   8.518538116671013902180962914473967738771E-1L,
 528   3.177729183645800174212539541058292579009E0L,
 529   4.006745668510308096259753538973038902990E0L,
 530   /* 1.000000000000000000000000000000000000000E0L, */
 531 };
 532
 533 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 534    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 535    2.85711669921875 <= x <= 4.54541015625
 536    Peak relative error 2.2e-21 */
 537 #ifdef __STDC__
 538 static const long double qr3[7] = {
 539 #else
 540 static long double qr3[7] = {
 541 #endif
 542   -3.618746299358445926506719188614570588404E-6L,
 543   -2.951146018465419674063882650970344502798E-4L,
 544   -7.728518171262562194043409753656506795258E-3L,
 545   -8.058010968753999435006488158237984014883E-2L,
 546   -3.356232856677966691703904770937143483472E-1L,
 547   -4.858192581793118040782557808823460276452E-1L,
 548   -1.592399251246473643510898335746432479373E-1L,
 549 };
 550 #ifdef __STDC__
 551 static const long double qs3[7] = {
 552 #else
 553 static long double qs3[7] = {
 554 #endif
 555   3.529139957987837084554591421329876744262E-5L,
 556   2.973602667215766676998703687065066180115E-3L,
 557   8.273534546240864308494062287908662592100E-2L,
 558   9.613359842126507198241321110649974032726E-1L,
 559   4.853923697093974370118387947065402707519E0L,
 560   1.002671608961669247462020977417828796933E1L,
 561   7.028927383922483728931327850683151410267E0L,
 562   /* 1.000000000000000000000000000000000000000E0L, */
 563 };
 564
 565 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 566    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 567    2 <= x <= 2.85711669921875
 568    Peak relative error 6.9e-22 */
 569 #ifdef __STDC__
 570 static const long double qr2[7] = {
 571 #else
 572 static long double qr2[7] = {
 573 #endif
 574   -1.372751603025230017220666013816502528318E-4L,
 575   -6.879190253347766576229143006767218972834E-3L,
 576   -1.061253572090925414598304855316280077828E-1L,
 577   -6.262164224345471241219408329354943337214E-1L,
 578   -1.423149636514768476376254324731437473915E0L,
 579   -1.087955310491078933531734062917489870754E0L,
 580   -1.826821119773182847861406108689273719137E-1L,
 581 };
 582 #ifdef __STDC__
 583 static const long double qs2[7] = {
 584 #else
 585 static long double qs2[7] = {
 586 #endif
 587   1.338768933634451601814048220627185324007E-3L,
 588   7.071099998918497559736318523932241901810E-2L,
 589   1.200511429784048632105295629933382142221E0L,
 590   8.327301713640367079030141077172031825276E0L,
 591   2.468479301872299311658145549931764426840E1L,
 592   2.961179686096262083509383820557051621644E1L,
 593   1.201402313144305153005639494661767354977E1L,
 594  /* 1.000000000000000000000000000000000000000E0L, */
 595 };
 596
 597
 598 #ifdef __STDC__
 599 static long double
 600 qone (long double x)
 601 #else
 602 static long double
 603 qone (x)
 604      long double x;
 605 #endif
 606 {
 607 #ifdef __STDC__
 608   const long double *p, *q;
 609 #else
 610   long double *p, *q;
 611 #endif
 612   static long double s, r, z;
 613   int32_t ix;
 614   u_int32_t se, i0, i1;
 615
 616   GET_LDOUBLE_WORDS (se, i0, i1, x);
 617   ix = se & 0x7fff;
 618   if (ix >= 0x4002)             /* x >= 8 */
 619     {
 620       p = qr8;
 621       q = qs8;
 622     }
 623   else
 624     {
 625       i1 = (ix << 16) | (i0 >> 16);
 626       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 627         {
 628           p = qr5;
 629           q = qs5;
 630         }
 631       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 632         {
 633           p = qr3;
 634           q = qs3;
 635         }
 636       else if (ix >= 0x4000)    /* x better be >= 2 */
 637         {
 638           p = qr2;
 639           q = qs2;
 640         }
 641     }
 642   z = one / (x * x);
 643   r =
 644     p[0] + z * (p[1] +
 645                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 646   s =
 647     q[0] + z * (q[1] +
 648                 z * (q[2] +
 649                      z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
 650   return (.375 + z * r / s) / x;
 651 }