Back to index

glibc  2.9
e_j1l.c
Go to the documentation of this file.
00001 /*
00002  * ====================================================
00003  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00004  *
00005  * Developed at SunPro, a Sun Microsystems, Inc. business.
00006  * Permission to use, copy, modify, and distribute this
00007  * software is freely granted, provided that this notice
00008  * is preserved.
00009  * ====================================================
00010  */
00011 
00012 /* Long double expansions are
00013   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
00014   and are incorporated herein by permission of the author.  The author 
00015   reserves the right to distribute this material elsewhere under different
00016   copying permissions.  These modifications are distributed here under 
00017   the following terms:
00018 
00019     This library is free software; you can redistribute it and/or
00020     modify it under the terms of the GNU Lesser General Public
00021     License as published by the Free Software Foundation; either
00022     version 2.1 of the License, or (at your option) any later version.
00023 
00024     This library is distributed in the hope that it will be useful,
00025     but WITHOUT ANY WARRANTY; without even the implied warranty of
00026     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00027     Lesser General Public License for more details.
00028 
00029     You should have received a copy of the GNU Lesser General Public
00030     License along with this library; if not, write to the Free Software
00031     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00032 
00033 /* __ieee754_j1(x), __ieee754_y1(x)
00034  * Bessel function of the first and second kinds of order zero.
00035  * Method -- j1(x):
00036  *     1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
00037  *     2. Reduce x to |x| since j1(x)=-j1(-x),  and
00038  *        for x in (0,2)
00039  *            j1(x) = x/2 + x*z*R0/S0,  where z = x*x;
00040  *        for x in (2,inf)
00041  *            j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
00042  *            y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
00043  *        where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
00044  *        as follow:
00045  *            cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
00046  *                   =  1/sqrt(2) * (sin(x) - cos(x))
00047  *            sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
00048  *                   = -1/sqrt(2) * (sin(x) + cos(x))
00049  *        (To avoid cancellation, use
00050  *            sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
00051  *         to compute the worse one.)
00052  *
00053  *     3 Special cases
00054  *            j1(nan)= nan
00055  *            j1(0) = 0
00056  *            j1(inf) = 0
00057  *
00058  * Method -- y1(x):
00059  *     1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
00060  *     2. For x<2.
00061  *        Since
00062  *            y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
00063  *        therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
00064  *        We use the following function to approximate y1,
00065  *            y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
00066  *        Note: For tiny x, 1/x dominate y1 and hence
00067  *            y1(tiny) = -2/pi/tiny
00068  *     3. For x>=2.
00069  *            y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
00070  *        where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
00071  *        by method mentioned above.
00072  */
00073 
00074 #include "math.h"
00075 #include "math_private.h"
00076 
00077 #ifdef __STDC__
00078 static long double pone (long double), qone (long double);
00079 #else
00080 static long double pone (), qone ();
00081 #endif
00082 
00083 #ifdef __STDC__
00084 static const long double
00085 #else
00086 static long double
00087 #endif
00088   huge = 1e4930L,
00089  one = 1.0L,
00090  invsqrtpi = 5.6418958354775628694807945156077258584405e-1L,
00091   tpi =  6.3661977236758134307553505349005744813784e-1L,
00092 
00093   /* J1(x) = .5 x + x x^2 R(x^2) / S(x^2)
00094      0 <= x <= 2
00095      Peak relative error 4.5e-21 */
00096 R[5] = {
00097     -9.647406112428107954753770469290757756814E7L,
00098     2.686288565865230690166454005558203955564E6L,
00099     -3.689682683905671185891885948692283776081E4L,
00100     2.195031194229176602851429567792676658146E2L,
00101     -5.124499848728030297902028238597308971319E-1L,
00102 },
00103 
00104   S[4] =
00105 {
00106   1.543584977988497274437410333029029035089E9L,
00107   2.133542369567701244002565983150952549520E7L,
00108   1.394077011298227346483732156167414670520E5L,
00109   5.252401789085732428842871556112108446506E2L,
00110   /*  1.000000000000000000000000000000000000000E0L, */
00111 };
00112 
00113 #ifdef __STDC__
00114 static const long double zero = 0.0;
00115 #else
00116 static long double zero = 0.0;
00117 #endif
00118 
00119 
00120 #ifdef __STDC__
00121 long double
00122 __ieee754_j1l (long double x)
00123 #else
00124 long double
00125 __ieee754_j1l (x)
00126      long double x;
00127 #endif
00128 {
00129   long double z, c, r, s, ss, cc, u, v, y;
00130   int32_t ix;
00131   u_int32_t se, i0, i1;
00132 
00133   GET_LDOUBLE_WORDS (se, i0, i1, x);
00134   ix = se & 0x7fff;
00135   if (ix >= 0x7fff)
00136     return one / x;
00137   y = fabsl (x);
00138   if (ix >= 0x4000)
00139     {                       /* |x| >= 2.0 */
00140       __sincosl (y, &s, &c);
00141       ss = -s - c;
00142       cc = s - c;
00143       if (ix < 0x7ffe)
00144        {                    /* make sure y+y not overflow */
00145          z = __cosl (y + y);
00146          if ((s * c) > zero)
00147            cc = z / ss;
00148          else
00149            ss = z / cc;
00150        }
00151       /*
00152        * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
00153        * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
00154        */
00155       if (ix > 0x4080)
00156        z = (invsqrtpi * cc) / __ieee754_sqrtl (y);
00157       else
00158        {
00159          u = pone (y);
00160          v = qone (y);
00161          z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (y);
00162        }
00163       if (se & 0x8000)
00164        return -z;
00165       else
00166        return z;
00167     }
00168   if (ix < 0x3fde) /* |x| < 2^-33 */
00169     {
00170       if (huge + x > one)
00171        return 0.5 * x;             /* inexact if x!=0 necessary */
00172     }
00173   z = x * x;
00174   r = z * (R[0] + z * (R[1]+ z * (R[2] + z * (R[3] + z * R[4]))));
00175   s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z)));
00176   r *= x;
00177   return (x * 0.5 + r / s);
00178 }
00179 
00180 
00181 /* Y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + x R(x^2)
00182    0 <= x <= 2
00183    Peak relative error 2.3e-23 */
00184 #ifdef __STDC__
00185 static const long double U0[6] = {
00186 #else
00187 static long double U0[6] = {
00188 #endif
00189   -5.908077186259914699178903164682444848615E10L,
00190   1.546219327181478013495975514375773435962E10L,
00191   -6.438303331169223128870035584107053228235E8L,
00192   9.708540045657182600665968063824819371216E6L,
00193   -6.138043997084355564619377183564196265471E4L,
00194   1.418503228220927321096904291501161800215E2L,
00195 };
00196 #ifdef __STDC__
00197 static const long double V0[5] = {
00198 #else
00199 static long double V0[5] = {
00200 #endif
00201   3.013447341682896694781964795373783679861E11L,
00202   4.669546565705981649470005402243136124523E9L,
00203   3.595056091631351184676890179233695857260E7L,
00204   1.761554028569108722903944659933744317994E5L,
00205   5.668480419646516568875555062047234534863E2L,
00206   /*  1.000000000000000000000000000000000000000E0L, */
00207 };
00208 
00209 
00210 #ifdef __STDC__
00211 long double
00212 __ieee754_y1l (long double x)
00213 #else
00214 long double
00215 __ieee754_y1l (x)
00216      long double x;
00217 #endif
00218 {
00219   long double z, s, c, ss, cc, u, v;
00220   int32_t ix;
00221   u_int32_t se, i0, i1;
00222 
00223   GET_LDOUBLE_WORDS (se, i0, i1, x);
00224   ix = se & 0x7fff;
00225   /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
00226   if (se & 0x8000)
00227     return zero / (zero * x);
00228   if (ix >= 0x7fff)
00229     return one / (x + x * x);
00230   if ((i0 | i1) == 0)
00231     return -HUGE_VALL + x;  /* -inf and overflow exception.  */
00232   if (ix >= 0x4000)
00233     {                       /* |x| >= 2.0 */
00234       __sincosl (x, &s, &c);
00235       ss = -s - c;
00236       cc = s - c;
00237       if (ix < 0x7fe00000)
00238        {                    /* make sure x+x not overflow */
00239          z = __cosl (x + x);
00240          if ((s * c) > zero)
00241            cc = z / ss;
00242          else
00243            ss = z / cc;
00244        }
00245       /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
00246        * where x0 = x-3pi/4
00247        *      Better formula:
00248        *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
00249        *                      =  1/sqrt(2) * (sin(x) - cos(x))
00250        *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
00251        *                      = -1/sqrt(2) * (cos(x) + sin(x))
00252        * To avoid cancellation, use
00253        *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
00254        * to compute the worse one.
00255        */
00256       if (ix > 0x4080)
00257        z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
00258       else
00259        {
00260          u = pone (x);
00261          v = qone (x);
00262          z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x);
00263        }
00264       return z;
00265     }
00266   if (ix <= 0x3fbe)
00267     {                       /* x < 2**-65 */
00268       return (-tpi / x);
00269     }
00270   z = x * x;
00271  u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * (U0[4] + z * U0[5]))));
00272  v = V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * (V0[4] + z))));
00273   return (x * (u / v) +
00274          tpi * (__ieee754_j1l (x) * __ieee754_logl (x) - one / x));
00275 }
00276 
00277 
00278 /* For x >= 8, the asymptotic expansions of pone is
00279  *     1 + 15/128 s^2 - 4725/2^15 s^4 - ...,     where s = 1/x.
00280  * We approximate pone by
00281  *     pone(x) = 1 + (R/S)
00282  */
00283 
00284 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
00285    P1(x) = 1 + z^2 R(z^2), z=1/x
00286    8 <= x <= inf  (0 <= z <= 0.125)
00287    Peak relative error 5.2e-22  */
00288 
00289 #ifdef __STDC__
00290 static const long double pr8[7] = {
00291 #else
00292 static long double pr8[7] = {
00293 #endif
00294   8.402048819032978959298664869941375143163E-9L,
00295   1.813743245316438056192649247507255996036E-6L,
00296   1.260704554112906152344932388588243836276E-4L,
00297   3.439294839869103014614229832700986965110E-3L,
00298   3.576910849712074184504430254290179501209E-2L,
00299   1.131111483254318243139953003461511308672E-1L,
00300   4.480715825681029711521286449131671880953E-2L,
00301 };
00302 #ifdef __STDC__
00303 static const long double ps8[6] = {
00304 #else
00305 static long double ps8[6] = {
00306 #endif
00307   7.169748325574809484893888315707824924354E-8L,
00308   1.556549720596672576431813934184403614817E-5L,
00309   1.094540125521337139209062035774174565882E-3L,
00310   3.060978962596642798560894375281428805840E-2L,
00311   3.374146536087205506032643098619414507024E-1L,
00312   1.253830208588979001991901126393231302559E0L,
00313   /* 1.000000000000000000000000000000000000000E0L, */
00314 };
00315 
00316 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
00317    P1(x) = 1 + z^2 R(z^2), z=1/x
00318    4.54541015625 <= x <= 8
00319    Peak relative error 7.7e-22  */
00320 #ifdef __STDC__
00321 static const long double pr5[7] = {
00322 #else
00323 static long double pr5[7] = {
00324 #endif
00325   4.318486887948814529950980396300969247900E-7L,
00326   4.715341880798817230333360497524173929315E-5L,
00327   1.642719430496086618401091544113220340094E-3L,
00328   2.228688005300803935928733750456396149104E-2L,
00329   1.142773760804150921573259605730018327162E-1L,
00330   1.755576530055079253910829652698703791957E-1L,
00331   3.218803858282095929559165965353784980613E-2L,
00332 };
00333 #ifdef __STDC__
00334 static const long double ps5[6] = {
00335 #else
00336 static long double ps5[6] = {
00337 #endif
00338   3.685108812227721334719884358034713967557E-6L,
00339   4.069102509511177498808856515005792027639E-4L,
00340   1.449728676496155025507893322405597039816E-2L,
00341   2.058869213229520086582695850441194363103E-1L,
00342   1.164890985918737148968424972072751066553E0L,
00343   2.274776933457009446573027260373361586841E0L,
00344   /*  1.000000000000000000000000000000000000000E0L,*/
00345 };
00346 
00347 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
00348    P1(x) = 1 + z^2 R(z^2), z=1/x
00349    2.85711669921875 <= x <= 4.54541015625
00350    Peak relative error 6.5e-21  */
00351 #ifdef __STDC__
00352 static const long double pr3[7] = {
00353 #else
00354 static long double pr3[7] = {
00355 #endif
00356   1.265251153957366716825382654273326407972E-5L,
00357   8.031057269201324914127680782288352574567E-4L,
00358   1.581648121115028333661412169396282881035E-2L,
00359   1.179534658087796321928362981518645033967E-1L,
00360   3.227936912780465219246440724502790727866E-1L,
00361   2.559223765418386621748404398017602935764E-1L,
00362   2.277136933287817911091370397134882441046E-2L,
00363 };
00364 #ifdef __STDC__
00365 static const long double ps3[6] = {
00366 #else
00367 static long double ps3[6] = {
00368 #endif
00369   1.079681071833391818661952793568345057548E-4L,
00370   6.986017817100477138417481463810841529026E-3L,
00371   1.429403701146942509913198539100230540503E-1L,
00372   1.148392024337075609460312658938700765074E0L,
00373   3.643663015091248720208251490291968840882E0L,
00374   3.990702269032018282145100741746633960737E0L,
00375   /* 1.000000000000000000000000000000000000000E0L, */
00376 };
00377 
00378 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
00379    P1(x) = 1 + z^2 R(z^2), z=1/x
00380    2 <= x <= 2.85711669921875
00381    Peak relative error 3.5e-21  */
00382 #ifdef __STDC__
00383 static const long double pr2[7] = {
00384 #else
00385 static long double pr2[7] = {
00386 #endif
00387   2.795623248568412225239401141338714516445E-4L,
00388   1.092578168441856711925254839815430061135E-2L,
00389   1.278024620468953761154963591853679640560E-1L,
00390   5.469680473691500673112904286228351988583E-1L,
00391   8.313769490922351300461498619045639016059E-1L,
00392   3.544176317308370086415403567097130611468E-1L,
00393   1.604142674802373041247957048801599740644E-2L,
00394 };
00395 #ifdef __STDC__
00396 static const long double ps2[6] = {
00397 #else
00398 static long double ps2[6] = {
00399 #endif
00400   2.385605161555183386205027000675875235980E-3L,
00401   9.616778294482695283928617708206967248579E-2L,
00402   1.195215570959693572089824415393951258510E0L,
00403   5.718412857897054829999458736064922974662E0L,
00404   1.065626298505499086386584642761602177568E1L,
00405   6.809140730053382188468983548092322151791E0L,
00406  /* 1.000000000000000000000000000000000000000E0L, */
00407 };
00408 
00409 
00410 #ifdef __STDC__
00411 static long double
00412 pone (long double x)
00413 #else
00414 static long double
00415 pone (x)
00416      long double x;
00417 #endif
00418 {
00419 #ifdef __STDC__
00420   const long double *p, *q;
00421 #else
00422   long double *p, *q;
00423 #endif
00424   long double z, r, s;
00425   int32_t ix;
00426   u_int32_t se, i0, i1;
00427 
00428   GET_LDOUBLE_WORDS (se, i0, i1, x);
00429   ix = se & 0x7fff;
00430   if (ix >= 0x4002) /* x >= 8 */
00431     {
00432       p = pr8;
00433       q = ps8;
00434     }
00435   else
00436     {
00437       i1 = (ix << 16) | (i0 >> 16);
00438       if (i1 >= 0x40019174) /* x >= 4.54541015625 */
00439        {
00440          p = pr5;
00441          q = ps5;
00442        }
00443       else if (i1 >= 0x4000b6db)   /* x >= 2.85711669921875 */
00444        {
00445          p = pr3;
00446          q = ps3;
00447        }
00448       else if (ix >= 0x4000)       /* x better be >= 2 */
00449        {
00450          p = pr2;
00451          q = ps2;
00452        }
00453     }
00454   z = one / (x * x);
00455  r = p[0] + z * (p[1] +
00456                z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
00457  s = q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
00458   return one + z * r / s;
00459 }
00460 
00461 
00462 /* For x >= 8, the asymptotic expansions of qone is
00463  *     3/8 s - 105/1024 s^3 - ..., where s = 1/x.
00464  * We approximate pone by
00465  *     qone(x) = s*(0.375 + (R/S))
00466  */
00467 
00468 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
00469    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
00470    8 <= x <= inf
00471    Peak relative error 8.3e-22 */
00472 
00473 #ifdef __STDC__
00474 static const long double qr8[7] = {
00475 #else
00476 static long double qr8[7] = {
00477 #endif
00478   -5.691925079044209246015366919809404457380E-10L,
00479   -1.632587664706999307871963065396218379137E-7L,
00480   -1.577424682764651970003637263552027114600E-5L,
00481   -6.377627959241053914770158336842725291713E-4L,
00482   -1.087408516779972735197277149494929568768E-2L,
00483   -6.854943629378084419631926076882330494217E-2L,
00484   -1.055448290469180032312893377152490183203E-1L,
00485 };
00486 #ifdef __STDC__
00487 static const long double qs8[7] = {
00488 #else
00489 static long double qs8[7] = {
00490 #endif
00491   5.550982172325019811119223916998393907513E-9L,
00492   1.607188366646736068460131091130644192244E-6L,
00493   1.580792530091386496626494138334505893599E-4L,
00494   6.617859900815747303032860443855006056595E-3L,
00495   1.212840547336984859952597488863037659161E-1L,
00496   9.017885953937234900458186716154005541075E-1L,
00497   2.201114489712243262000939120146436167178E0L,
00498   /* 1.000000000000000000000000000000000000000E0L, */
00499 };
00500 
00501 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
00502    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
00503    4.54541015625 <= x <= 8
00504    Peak relative error 4.1e-22 */
00505 #ifdef __STDC__
00506 static const long double qr5[7] = {
00507 #else
00508 static long double qr5[7] = {
00509 #endif
00510   -6.719134139179190546324213696633564965983E-8L,
00511   -9.467871458774950479909851595678622044140E-6L,
00512   -4.429341875348286176950914275723051452838E-4L,
00513   -8.539898021757342531563866270278505014487E-3L,
00514   -6.818691805848737010422337101409276287170E-2L,
00515   -1.964432669771684034858848142418228214855E-1L,
00516   -1.333896496989238600119596538299938520726E-1L,
00517 };
00518 #ifdef __STDC__
00519 static const long double qs5[7] = {
00520 #else
00521 static long double qs5[7] = {
00522 #endif
00523   6.552755584474634766937589285426911075101E-7L,
00524   9.410814032118155978663509073200494000589E-5L,
00525   4.561677087286518359461609153655021253238E-3L,
00526   9.397742096177905170800336715661091535805E-2L,
00527   8.518538116671013902180962914473967738771E-1L,
00528   3.177729183645800174212539541058292579009E0L,
00529   4.006745668510308096259753538973038902990E0L,
00530   /* 1.000000000000000000000000000000000000000E0L, */
00531 };
00532 
00533 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
00534    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
00535    2.85711669921875 <= x <= 4.54541015625
00536    Peak relative error 2.2e-21 */
00537 #ifdef __STDC__
00538 static const long double qr3[7] = {
00539 #else
00540 static long double qr3[7] = {
00541 #endif
00542   -3.618746299358445926506719188614570588404E-6L,
00543   -2.951146018465419674063882650970344502798E-4L,
00544   -7.728518171262562194043409753656506795258E-3L,
00545   -8.058010968753999435006488158237984014883E-2L,
00546   -3.356232856677966691703904770937143483472E-1L,
00547   -4.858192581793118040782557808823460276452E-1L,
00548   -1.592399251246473643510898335746432479373E-1L,
00549 };
00550 #ifdef __STDC__
00551 static const long double qs3[7] = {
00552 #else
00553 static long double qs3[7] = {
00554 #endif
00555   3.529139957987837084554591421329876744262E-5L,
00556   2.973602667215766676998703687065066180115E-3L,
00557   8.273534546240864308494062287908662592100E-2L,
00558   9.613359842126507198241321110649974032726E-1L,
00559   4.853923697093974370118387947065402707519E0L,
00560   1.002671608961669247462020977417828796933E1L,
00561   7.028927383922483728931327850683151410267E0L,
00562   /* 1.000000000000000000000000000000000000000E0L, */
00563 };
00564 
00565 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
00566    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
00567    2 <= x <= 2.85711669921875
00568    Peak relative error 6.9e-22 */
00569 #ifdef __STDC__
00570 static const long double qr2[7] = {
00571 #else
00572 static long double qr2[7] = {
00573 #endif
00574   -1.372751603025230017220666013816502528318E-4L,
00575   -6.879190253347766576229143006767218972834E-3L,
00576   -1.061253572090925414598304855316280077828E-1L,
00577   -6.262164224345471241219408329354943337214E-1L,
00578   -1.423149636514768476376254324731437473915E0L,
00579   -1.087955310491078933531734062917489870754E0L,
00580   -1.826821119773182847861406108689273719137E-1L,
00581 };
00582 #ifdef __STDC__
00583 static const long double qs2[7] = {
00584 #else
00585 static long double qs2[7] = {
00586 #endif
00587   1.338768933634451601814048220627185324007E-3L,
00588   7.071099998918497559736318523932241901810E-2L,
00589   1.200511429784048632105295629933382142221E0L,
00590   8.327301713640367079030141077172031825276E0L,
00591   2.468479301872299311658145549931764426840E1L,
00592   2.961179686096262083509383820557051621644E1L,
00593   1.201402313144305153005639494661767354977E1L,
00594  /* 1.000000000000000000000000000000000000000E0L, */
00595 };
00596 
00597 
00598 #ifdef __STDC__
00599 static long double
00600 qone (long double x)
00601 #else
00602 static long double
00603 qone (x)
00604      long double x;
00605 #endif
00606 {
00607 #ifdef __STDC__
00608   const long double *p, *q;
00609 #else
00610   long double *p, *q;
00611 #endif
00612   static long double s, r, z;
00613   int32_t ix;
00614   u_int32_t se, i0, i1;
00615 
00616   GET_LDOUBLE_WORDS (se, i0, i1, x);
00617   ix = se & 0x7fff;
00618   if (ix >= 0x4002)         /* x >= 8 */
00619     {
00620       p = qr8;
00621       q = qs8;
00622     }
00623   else
00624     {
00625       i1 = (ix << 16) | (i0 >> 16);
00626       if (i1 >= 0x40019174) /* x >= 4.54541015625 */
00627        {
00628          p = qr5;
00629          q = qs5;
00630        }
00631       else if (i1 >= 0x4000b6db)   /* x >= 2.85711669921875 */
00632        {
00633          p = qr3;
00634          q = qs3;
00635        }
00636       else if (ix >= 0x4000)       /* x better be >= 2 */
00637        {
00638          p = qr2;
00639          q = qs2;
00640        }
00641     }
00642   z = one / (x * x);
00643   r =
00644     p[0] + z * (p[1] +
00645               z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
00646   s =
00647     q[0] + z * (q[1] +
00648               z * (q[2] +
00649                    z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
00650   return (.375 + z * r / s) / x;
00651 }