Back to index

glibc  2.9
e_j0l.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_j0(x), __ieee754_y0(x)
00034  * Bessel function of the first and second kinds of order zero.
00035  * Method -- j0(x):
00036  *     1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
00037  *     2. Reduce x to |x| since j0(x)=j0(-x),  and
00038  *        for x in (0,2)
00039  *            j0(x) = 1 - z/4 + z^2*R0/S0,  where z = x*x;
00040  *        for x in (2,inf)
00041  *            j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
00042  *        where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
00043  *        as follow:
00044  *            cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
00045  *                   = 1/sqrt(2) * (cos(x) + sin(x))
00046  *            sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
00047  *                   = 1/sqrt(2) * (sin(x) - cos(x))
00048  *        (To avoid cancellation, use
00049  *            sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
00050  *         to compute the worse one.)
00051  *
00052  *     3 Special cases
00053  *            j0(nan)= nan
00054  *            j0(0) = 1
00055  *            j0(inf) = 0
00056  *
00057  * Method -- y0(x):
00058  *     1. For x<2.
00059  *        Since
00060  *            y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
00061  *        therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
00062  *        We use the following function to approximate y0,
00063  *            y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
00064  *
00065  *        Note: For tiny x, U/V = u0 and j0(x)~1, hence
00066  *            y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
00067  *     2. For x>=2.
00068  *            y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
00069  *        where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
00070  *        by the method mentioned above.
00071  *     3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
00072  */
00073 
00074 #include "math.h"
00075 #include "math_private.h"
00076 
00077 #ifdef __STDC__
00078 static long double pzero (long double), qzero (long double);
00079 #else
00080 static long double pzero (), qzero ();
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   /* J0(x) = 1 - x^2 / 4 + x^4 R0(x^2) / S0(x^2)
00094      0 <= x <= 2
00095      peak relative error 1.41e-22 */
00096   R[5] = {
00097   4.287176872744686992880841716723478740566E7L,
00098   -6.652058897474241627570911531740907185772E5L,
00099   7.011848381719789863458364584613651091175E3L,
00100   -3.168040850193372408702135490809516253693E1L,
00101   6.030778552661102450545394348845599300939E-2L,
00102 },
00103 
00104  S[4] = {
00105    2.743793198556599677955266341699130654342E9L,
00106    3.364330079384816249840086842058954076201E7L,
00107    1.924119649412510777584684927494642526573E5L,
00108    6.239282256012734914211715620088714856494E2L,
00109    /*   1.000000000000000000000000000000000000000E0L,*/
00110 };
00111 
00112 #ifdef __STDC__
00113 static const long double zero = 0.0;
00114 #else
00115 static long double zero = 0.0;
00116 #endif
00117 
00118 #ifdef __STDC__
00119 long double
00120 __ieee754_j0l (long double x)
00121 #else
00122 long double
00123 __ieee754_j0l (x)
00124      long double x;
00125 #endif
00126 {
00127   long double z, s, c, ss, cc, r, u, v;
00128   int32_t ix;
00129   u_int32_t se, i0, i1;
00130 
00131   GET_LDOUBLE_WORDS (se, i0, i1, x);
00132   ix = se & 0x7fff;
00133   if (ix >= 0x7fff)
00134     return one / (x * x);
00135   x = fabsl (x);
00136   if (ix >= 0x4000)         /* |x| >= 2.0 */
00137     {
00138       __sincosl (x, &s, &c);
00139       ss = s - c;
00140       cc = s + c;
00141       if (ix < 0x7ffe)
00142        {                    /* make sure x+x not overflow */
00143          z = -__cosl (x + x);
00144          if ((s * c) < zero)
00145            cc = z / ss;
00146          else
00147            ss = z / cc;
00148        }
00149       /*
00150        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
00151        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
00152        */
00153       if (ix > 0x4080)      /* 2^129 */
00154        z = (invsqrtpi * cc) / __ieee754_sqrtl (x);
00155       else
00156        {
00157          u = pzero (x);
00158          v = qzero (x);
00159          z = invsqrtpi * (u * cc - v * ss) / __ieee754_sqrtl (x);
00160        }
00161       return z;
00162     }
00163   if (ix < 0x3fef) /* |x| < 2**-16 */
00164     {
00165       if (huge + x > one)
00166        {                    /* raise inexact if x != 0 */
00167          if (ix < 0x3fde) /* |x| < 2^-33 */
00168            return one;
00169          else
00170            return one - 0.25 * x * x;
00171        }
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   if (ix < 0x3fff)
00177     {                       /* |x| < 1.00 */
00178       return (one - 0.25 * z + z * (r / s));
00179     }
00180   else
00181     {
00182       u = 0.5 * x;
00183       return ((one + u) * (one - u) + z * (r / s));
00184     }
00185 }
00186 
00187 
00188 /* y0(x) = 2/pi ln(x) J0(x) + U(x^2)/V(x^2)
00189    0 < x <= 2
00190    peak relative error 1.7e-21 */
00191 #ifdef __STDC__
00192 static const long double
00193 #else
00194 static long double
00195 #endif
00196 U[6] = {
00197   -1.054912306975785573710813351985351350861E10L,
00198   2.520192609749295139432773849576523636127E10L,
00199   -1.856426071075602001239955451329519093395E9L,
00200   4.079209129698891442683267466276785956784E7L,
00201   -3.440684087134286610316661166492641011539E5L,
00202   1.005524356159130626192144663414848383774E3L,
00203 };
00204 #ifdef __STDC__
00205 static const long double
00206 #else
00207 static long double
00208 #endif
00209 V[5] = {
00210   1.429337283720789610137291929228082613676E11L,
00211   2.492593075325119157558811370165695013002E9L,
00212   2.186077620785925464237324417623665138376E7L,
00213   1.238407896366385175196515057064384929222E5L,
00214   4.693924035211032457494368947123233101664E2L,
00215   /*  1.000000000000000000000000000000000000000E0L */
00216 };
00217 
00218 #ifdef __STDC__
00219 long double
00220 __ieee754_y0l (long double x)
00221 #else
00222 long double
00223 __ieee754_y0l (x)
00224      long double x;
00225 #endif
00226 {
00227   long double z, s, c, ss, cc, u, v;
00228   int32_t ix;
00229   u_int32_t se, i0, i1;
00230 
00231   GET_LDOUBLE_WORDS (se, i0, i1, x);
00232   ix = se & 0x7fff;
00233   /* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0  */
00234   if (se & 0x8000)
00235     return zero / (zero * x);
00236   if (ix >= 0x7fff)
00237     return one / (x + x * x);
00238   if ((i0 | i1) == 0)
00239     return -HUGE_VALL + x;  /* -inf and overflow exception.  */
00240   if (ix >= 0x4000)
00241     {                       /* |x| >= 2.0 */
00242 
00243       /* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
00244        * where x0 = x-pi/4
00245        *      Better formula:
00246        *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
00247        *                      =  1/sqrt(2) * (sin(x) + cos(x))
00248        *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
00249        *                      =  1/sqrt(2) * (sin(x) - cos(x))
00250        * To avoid cancellation, use
00251        *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
00252        * to compute the worse one.
00253        */
00254       __sincosl (x, &s, &c);
00255       ss = s - c;
00256       cc = s + c;
00257       /*
00258        * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
00259        * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
00260        */
00261       if (ix < 0x7ffe)
00262        {                    /* make sure x+x not overflow */
00263          z = -__cosl (x + x);
00264          if ((s * c) < zero)
00265            cc = z / ss;
00266          else
00267            ss = z / cc;
00268        }
00269       if (ix > 0x4080)      /* 1e39 */
00270        z = (invsqrtpi * ss) / __ieee754_sqrtl (x);
00271       else
00272        {
00273          u = pzero (x);
00274          v = qzero (x);
00275          z = invsqrtpi * (u * ss + v * cc) / __ieee754_sqrtl (x);
00276        }
00277       return z;
00278     }
00279   if (ix <= 0x3fde) /* x < 2^-33 */
00280     {
00281       z = -7.380429510868722527629822444004602747322E-2L
00282        + tpi * __ieee754_logl (x);
00283       return z;
00284     }
00285   z = x * x;
00286   u = U[0] + z * (U[1] + z * (U[2] + z * (U[3] + z * (U[4] + z * U[5]))));
00287   v = V[0] + z * (V[1] + z * (V[2] + z * (V[3] + z * (V[4] + z))));
00288   return (u / v + tpi * (__ieee754_j0l (x) * __ieee754_logl (x)));
00289 }
00290 
00291 /* The asymptotic expansions of pzero is
00292  *     1 - 9/128 s^2 + 11025/98304 s^4 - ...,    where s = 1/x.
00293  * For x >= 2, We approximate pzero by
00294  *     pzero(x) = 1 + s^2 R(s^2) / S(s^2)
00295  */
00296 #ifdef __STDC__
00297 static const long double pR8[7] = {
00298 #else
00299 static long double pR8[7] = {
00300 #endif
00301   /* 8 <= x <= inf
00302      Peak relative error 4.62 */
00303   -4.094398895124198016684337960227780260127E-9L,
00304   -8.929643669432412640061946338524096893089E-7L,
00305   -6.281267456906136703868258380673108109256E-5L,
00306   -1.736902783620362966354814353559382399665E-3L,
00307   -1.831506216290984960532230842266070146847E-2L,
00308   -5.827178869301452892963280214772398135283E-2L,
00309   -2.087563267939546435460286895807046616992E-2L,
00310 };
00311 #ifdef __STDC__
00312 static const long double pS8[6] = {
00313 #else
00314 static long double pS8[6] = {
00315 #endif
00316   5.823145095287749230197031108839653988393E-8L,
00317   1.279281986035060320477759999428992730280E-5L,
00318   9.132668954726626677174825517150228961304E-4L,
00319   2.606019379433060585351880541545146252534E-2L,
00320   2.956262215119520464228467583516287175244E-1L,
00321   1.149498145388256448535563278632697465675E0L,
00322   /* 1.000000000000000000000000000000000000000E0L, */
00323 };
00324 
00325 #ifdef __STDC__
00326 static const long double pR5[7] = {
00327 #else
00328 static long double pR5[7] = {
00329 #endif
00330   /* 4.54541015625 <= x <= 8
00331      Peak relative error 6.51E-22 */
00332   -2.041226787870240954326915847282179737987E-7L,
00333   -2.255373879859413325570636768224534428156E-5L,
00334   -7.957485746440825353553537274569102059990E-4L,
00335   -1.093205102486816696940149222095559439425E-2L,
00336   -5.657957849316537477657603125260701114646E-2L,
00337   -8.641175552716402616180994954177818461588E-2L,
00338   -1.354654710097134007437166939230619726157E-2L,
00339 };
00340 #ifdef __STDC__
00341 static const long double pS5[6] = {
00342 #else
00343 static long double pS5[6] = {
00344 #endif
00345   2.903078099681108697057258628212823545290E-6L,
00346   3.253948449946735405975737677123673867321E-4L,
00347   1.181269751723085006534147920481582279979E-2L,
00348   1.719212057790143888884745200257619469363E-1L,
00349   1.006306498779212467670654535430694221924E0L,
00350   2.069568808688074324555596301126375951502E0L,
00351   /* 1.000000000000000000000000000000000000000E0L, */
00352 };
00353 
00354 #ifdef __STDC__
00355 static const long double pR3[7] = {
00356 #else
00357 static long double pR3[7] = {
00358 #endif
00359   /* 2.85711669921875 <= x <= 4.54541015625
00360      peak relative error 5.25e-21 */
00361   -5.755732156848468345557663552240816066802E-6L,
00362   -3.703675625855715998827966962258113034767E-4L,
00363   -7.390893350679637611641350096842846433236E-3L,
00364   -5.571922144490038765024591058478043873253E-2L,
00365   -1.531290690378157869291151002472627396088E-1L,
00366   -1.193350853469302941921647487062620011042E-1L,
00367   -8.567802507331578894302991505331963782905E-3L,
00368 };
00369 #ifdef __STDC__
00370 static const long double pS3[6] = {
00371 #else
00372 static long double pS3[6] = {
00373 #endif
00374   8.185931139070086158103309281525036712419E-5L,
00375   5.398016943778891093520574483111255476787E-3L,
00376   1.130589193590489566669164765853409621081E-1L,
00377   9.358652328786413274673192987670237145071E-1L,
00378   3.091711512598349056276917907005098085273E0L,
00379   3.594602474737921977972586821673124231111E0L,
00380   /* 1.000000000000000000000000000000000000000E0L, */
00381 };
00382 
00383 #ifdef __STDC__
00384 static const long double pR2[7] = {
00385 #else
00386 static long double pR2[7] = {
00387 #endif
00388   /* 2 <= x <= 2.85711669921875
00389      peak relative error 2.64e-21 */
00390   -1.219525235804532014243621104365384992623E-4L,
00391   -4.838597135805578919601088680065298763049E-3L,
00392   -5.732223181683569266223306197751407418301E-2L,
00393   -2.472947430526425064982909699406646503758E-1L,
00394   -3.753373645974077960207588073975976327695E-1L,
00395   -1.556241316844728872406672349347137975495E-1L,
00396   -5.355423239526452209595316733635519506958E-3L,
00397 };
00398 #ifdef __STDC__
00399 static const long double pS2[6] = {
00400 #else
00401 static long double pS2[6] = {
00402 #endif
00403   1.734442793664291412489066256138894953823E-3L,
00404   7.158111826468626405416300895617986926008E-2L,
00405   9.153839713992138340197264669867993552641E-1L,
00406   4.539209519433011393525841956702487797582E0L,
00407   8.868932430625331650266067101752626253644E0L,
00408   6.067161890196324146320763844772857713502E0L,
00409   /* 1.000000000000000000000000000000000000000E0L, */
00410 };
00411 
00412 #ifdef __STDC__
00413 static long double
00414 pzero (long double x)
00415 #else
00416 static long double
00417 pzero (x)
00418      long double x;
00419 #endif
00420 {
00421 #ifdef __STDC__
00422   const long double *p, *q;
00423 #else
00424   long double *p, *q;
00425 #endif
00426   long double z, r, s;
00427   int32_t ix;
00428   u_int32_t se, i0, i1;
00429 
00430   GET_LDOUBLE_WORDS (se, i0, i1, x);
00431   ix = se & 0x7fff;
00432   if (ix >= 0x4002)
00433     {
00434       p = pR8;
00435       q = pS8;
00436     }                       /* x >= 8 */
00437   else
00438     {
00439       i1 = (ix << 16) | (i0 >> 16);
00440       if (i1 >= 0x40019174) /* x >= 4.54541015625 */
00441        {
00442          p = pR5;
00443          q = pS5;
00444        }
00445       else if (i1 >= 0x4000b6db)   /* x >= 2.85711669921875 */
00446        {
00447          p = pR3;
00448          q = pS3;
00449        }
00450       else if (ix >= 0x4000)       /* x better be >= 2 */
00451        {
00452          p = pR2;
00453          q = pS2;
00454        }
00455     }
00456   z = one / (x * x);
00457   r =
00458     p[0] + z * (p[1] +
00459               z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
00460   s =
00461     q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
00462   return (one + z * r / s);
00463 }
00464 
00465 
00466 /* For x >= 8, the asymptotic expansions of qzero is
00467  *     -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
00468  * We approximate qzero by
00469  *     qzero(x) = s*(-.125 + R(s^2) / S(s^2))
00470  */
00471 #ifdef __STDC__
00472 static const long double qR8[7] = {
00473 #else
00474 static long double qR8[7] = {
00475 #endif
00476   /* 8 <= x <= inf
00477      peak relative error 2.23e-21 */
00478   3.001267180483191397885272640777189348008E-10L,
00479   8.693186311430836495238494289942413810121E-8L,
00480   8.496875536711266039522937037850596580686E-6L,
00481   3.482702869915288984296602449543513958409E-4L,
00482   6.036378380706107692863811938221290851352E-3L,
00483   3.881970028476167836382607922840452192636E-2L,
00484   6.132191514516237371140841765561219149638E-2L,
00485 };
00486 #ifdef __STDC__
00487 static const long double qS8[7] = {
00488 #else
00489 static long double qS8[7] = {
00490 #endif
00491   4.097730123753051126914971174076227600212E-9L,
00492   1.199615869122646109596153392152131139306E-6L,
00493   1.196337580514532207793107149088168946451E-4L,
00494   5.099074440112045094341500497767181211104E-3L,
00495   9.577420799632372483249761659674764460583E-2L,
00496   7.385243015344292267061953461563695918646E-1L,
00497   1.917266424391428937962682301561699055943E0L,
00498   /* 1.000000000000000000000000000000000000000E0L, */
00499 };
00500 
00501 #ifdef __STDC__
00502 static const long double qR5[7] = {
00503 #else
00504 static long double qR5[7] = {
00505 #endif
00506   /* 4.54541015625 <= x <= 8
00507      peak relative error 1.03e-21 */
00508   3.406256556438974327309660241748106352137E-8L,
00509   4.855492710552705436943630087976121021980E-6L,
00510   2.301011739663737780613356017352912281980E-4L,
00511   4.500470249273129953870234803596619899226E-3L,
00512   3.651376459725695502726921248173637054828E-2L,
00513   1.071578819056574524416060138514508609805E-1L,
00514   7.458950172851611673015774675225656063757E-2L,
00515 };
00516 #ifdef __STDC__
00517 static const long double qS5[7] = {
00518 #else
00519 static long double qS5[7] = {
00520 #endif
00521   4.650675622764245276538207123618745150785E-7L,
00522   6.773573292521412265840260065635377164455E-5L,
00523   3.340711249876192721980146877577806687714E-3L,
00524   7.036218046856839214741678375536970613501E-2L,
00525   6.569599559163872573895171876511377891143E-1L,
00526   2.557525022583599204591036677199171155186E0L,
00527   3.457237396120935674982927714210361269133E0L,
00528   /* 1.000000000000000000000000000000000000000E0L,*/
00529 };
00530 
00531 #ifdef __STDC__
00532 static const long double qR3[7] = {
00533 #else
00534 static long double qR3[7] = {
00535 #endif
00536   /* 2.85711669921875 <= x <= 4.54541015625
00537      peak relative error 5.24e-21 */
00538   1.749459596550816915639829017724249805242E-6L,
00539   1.446252487543383683621692672078376929437E-4L,
00540   3.842084087362410664036704812125005761859E-3L,
00541   4.066369994699462547896426554180954233581E-2L,
00542   1.721093619117980251295234795188992722447E-1L,
00543   2.538595333972857367655146949093055405072E-1L,
00544   8.560591367256769038905328596020118877936E-2L,
00545 };
00546 #ifdef __STDC__
00547 static const long double qS3[7] = {
00548 #else
00549 static long double qS3[7] = {
00550 #endif
00551   2.388596091707517488372313710647510488042E-5L,
00552   2.048679968058758616370095132104333998147E-3L,
00553   5.824663198201417760864458765259945181513E-2L,
00554   6.953906394693328750931617748038994763958E-1L,
00555   3.638186936390881159685868764832961092476E0L,
00556   7.900169524705757837298990558459547842607E0L,
00557   5.992718532451026507552820701127504582907E0L,
00558   /* 1.000000000000000000000000000000000000000E0L, */
00559 };
00560 
00561 #ifdef __STDC__
00562 static const long double qR2[7] = {
00563 #else
00564 static long double qR2[7] = {
00565 #endif
00566   /* 2 <= x <= 2.85711669921875
00567      peak relative error 1.58e-21  */
00568   6.306524405520048545426928892276696949540E-5L,
00569   3.209606155709930950935893996591576624054E-3L,
00570   5.027828775702022732912321378866797059604E-2L,
00571   3.012705561838718956481911477587757845163E-1L,
00572   6.960544893905752937420734884995688523815E-1L,
00573   5.431871999743531634887107835372232030655E-1L,
00574   9.447736151202905471899259026430157211949E-2L,
00575 };
00576 #ifdef __STDC__
00577 static const long double qS2[7] = {
00578 #else
00579 static long double qS2[7] = {
00580 #endif
00581   8.610579901936193494609755345106129102676E-4L,
00582   4.649054352710496997203474853066665869047E-2L,
00583   8.104282924459837407218042945106320388339E-1L,
00584   5.807730930825886427048038146088828206852E0L,
00585   1.795310145936848873627710102199881642939E1L,
00586   2.281313316875375733663657188888110605044E1L,
00587   1.011242067883822301487154844458322200143E1L,
00588   /* 1.000000000000000000000000000000000000000E0L, */
00589 };
00590 
00591 #ifdef __STDC__
00592 static long double
00593 qzero (long double x)
00594 #else
00595 static long double
00596 qzero (x)
00597      long double x;
00598 #endif
00599 {
00600 #ifdef __STDC__
00601   const long double *p, *q;
00602 #else
00603   long double *p, *q;
00604 #endif
00605   long double s, r, z;
00606   int32_t ix;
00607   u_int32_t se, i0, i1;
00608 
00609   GET_LDOUBLE_WORDS (se, i0, i1, x);
00610   ix = se & 0x7fff;
00611   if (ix >= 0x4002)         /* x >= 8 */
00612     {
00613       p = qR8;
00614       q = qS8;
00615     }
00616   else
00617     {
00618       i1 = (ix << 16) | (i0 >> 16);
00619       if (i1 >= 0x40019174) /* x >= 4.54541015625 */
00620        {
00621          p = qR5;
00622          q = qS5;
00623        }
00624       else if (i1 >= 0x4000b6db)   /* x >= 2.85711669921875 */
00625        {
00626          p = qR3;
00627          q = qS3;
00628        }
00629       else if (ix >= 0x4000)       /* x better be >= 2 */
00630        {
00631          p = qR2;
00632          q = qS2;
00633        }
00634     }
00635   z = one / (x * x);
00636   r =
00637     p[0] + z * (p[1] +
00638               z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
00639   s =
00640     q[0] + z * (q[1] +
00641               z * (q[2] +
00642                    z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
00643   return (-.125 + z * r / s) / x;
00644 }