Print this page
11210 libm should be cstyle(1ONBLD) clean


   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */

  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29 
  30 /*
  31  * Floating point Bessel's function of the first and second kinds
  32  * of order zero: j0(x),y0(x);
  33  *
  34  * Special cases:
  35  *      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  36  *      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
  37  */
  38 
  39 #pragma weak __j0l = j0l
  40 #pragma weak __y0l = y0l
  41 
  42 #include "libm.h"
  43 #include "longdouble.h"
  44 
  45 #define GENERIC long double
  46 static const GENERIC
  47 zero    = 0.0L,
  48 small   = 1.0e-9L,
  49 tiny    = 1.0e-38L,
  50 one     = 1.0L,
  51 five    = 5.0L,
  52 eight   = 8.0L,
  53 invsqrtpi= 5.641895835477562869480794515607725858441e-0001L,
  54 tpi     = 0.636619772367581343075535053490057448L;
  55 
  56 static GENERIC pzero(GENERIC);
  57 static GENERIC qzero(GENERIC);
  58 
  59 static GENERIC r0[7] = {
  60   -2.499999999999999999999999999999998934492e-0001L,
  61    1.272657927360049786327618451133763714880e-0002L,
  62   -2.694499763712963276900636693400659600898e-0004L,
  63    2.724877475058977576903234070919616447883e-0006L,
  64   -1.432617103214330236967477495393076320281e-0008L,
  65    3.823248804080079168706683540513792224471e-0011L,
  66   -4.183174277567983647337568504286313665065e-0014L,
  67 };

  68 static GENERIC s0[7] = {
  69    1.0e0L,
  70    1.159368290559800854689526195462884666395e-0002L,
  71    6.629397597394973383009743876169946772559e-0005L,
  72    2.426779981394054406305431142501735094340e-0007L,
  73    6.097663491248511069094400469635449749883e-0010L,
  74    1.017019133340929220238747413216052224036e-0012L,
  75    9.012593179306197579518374581969371278481e-0016L,
  76 };
  77 
  78 GENERIC
  79 j0l(x) GENERIC x;{
  80         GENERIC z, s,c,ss,cc,r,u,v;

  81         int i;
  82 
  83         if (isnanl(x)) return x+x;


  84         x = fabsl(x);

  85         if (x > 1.28L) {
  86                 if (!finitel(x)) return zero;


  87                 s = sinl(x);
  88                 c = cosl(x);


  89         /*
  90          * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
  91          * where x0 = x-pi/4
  92          *      Better formula:
  93          *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
  94          *                      = 1/sqrt(2) * (cos(x) + sin(x))
  95          *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
  96          *                      = 1/sqrt(2) * (sin(x) - cos(x))
  97          * To avoid cancellation, use
  98          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  99          * to compute the worse one.
 100          */
 101                 if (x>1.0e2450L) {   /* x+x may overflow */
 102                         ss = s-c;
 103                         cc = s+c;

 104                 } else if (signbitl(s) != signbitl(c)) {
 105                         ss = s - c;
 106                         cc = -cosl(x+x)/ss;
 107                 } else {
 108                         cc = s + c;
 109                         ss = -cosl(x+x)/cc;
 110                 }

 111         /*
 112          * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 113          * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 114          */
 115                 if (x>1.0e120L) return (invsqrtpi*cc)/sqrtl(x);
 116                 u = pzero(x); v = qzero(x);
 117                 return invsqrtpi*(u*cc-v*ss)/sqrtl(x);
 118         }
 119         if (x<=small) {
 120             if (x<=tiny) return one-x;
 121             else return one-x*x*0.25L;
 122         }
 123         z = x*x;
 124         r = r0[6]; s = s0[6];
 125         for(i=5;i>=0;i--) {
 126             r = r*z + r0[i];
 127             s = s*z + s0[i];
 128         }
 129         return(one+z*(r/s));










 130 }
 131 
 132 static const GENERIC u0[8] = {
 133   -7.380429510868722527434392794848301631220e-0002L,
 134    1.766855559625940791857536949301981816513e-0001L,
 135   -1.386470722701047923235553251240162839408e-0002L,
 136    3.520149242724811578636970811631224862615e-0004L,
 137   -3.978599663243790049853642275624951870025e-0006L,
 138    2.228801153263957224547222556806915479763e-0008L,
 139   -6.121246764298785018658597179498837316177e-0011L,
 140    6.677103629722678833475965810525587396596e-0014L,
 141 };

 142 static const GENERIC v0[8] = {
 143    1.0e0L,
 144    1.247164416539111311571676766127767127970e-0002L,
 145    7.829144749639791500052900281489367443576e-0005L,
 146    3.247126540422245330511218321013360336606e-0007L,
 147    9.750516724789499678567062572549568447869e-0010L,
 148    2.156713223173591212250543390258458098776e-0012L,
 149    3.322169561597890004231482431236452752624e-0015L,
 150    2.821213295314000924252226486305726805093e-0018L,
 151 };
 152 
 153 GENERIC
 154 y0l(x) GENERIC x;{
 155         GENERIC z, s,c,ss,cc,u,v;

 156         int i;
 157         volatile GENERIC d;
 158 
 159         if (isnanl(x)) return x+x;


 160         if (x <= zero) {
 161                 if (x == zero)
 162                     d= -one/(x-x);
 163                 else
 164                     d = zero/(x-x);
 165         }

 166 #ifdef lint
 167         d = d;
 168 #endif

 169         if (x > 1.28L) {
 170                 if (!finitel(x)) return zero;


 171                 s = sinl(x);
 172                 c = cosl(x);


 173         /*
 174          * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
 175          * where x0 = x-pi/4
 176          *      Better formula:
 177          *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
 178          *                      = 1/sqrt(2) * (cos(x) + sin(x))
 179          *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
 180          *                      = 1/sqrt(2) * (sin(x) - cos(x))
 181          * To avoid cancellation, use
 182          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 183          * to compute the worse one.
 184          */
 185                 if (x>1.0e2450L) {   /* x+x may overflow */
 186                         ss = s-c;
 187                         cc = s+c;

 188                 } else if (signbitl(s) != signbitl(c)) {
 189                         ss = s - c;
 190                         cc = -cosl(x+x)/ss;
 191                 } else {
 192                         cc = s + c;
 193                         ss = -cosl(x+x)/cc;
 194                 }

 195         /*
 196          * j0(x) = 1/sqrt(pi*x) * (P(0,x)*cc - Q(0,x)*ss)
 197          * y0(x) = 1/sqrt(pi*x) * (P(0,x)*ss + Q(0,x)*cc)
 198          */
 199                 if (x>1.0e120L) return (invsqrtpi*ss)/sqrtl(x);
 200                 return invsqrtpi*(pzero(x)*ss+qzero(x)*cc)/sqrtl(x);
 201 

 202         }
 203         if (x<=tiny) {
 204             return (u0[0] + tpi*logl(x));
 205         }
 206         z = x*x;
 207         u = u0[7]; v = v0[7];
 208         for(i=6;i>=0;i--) {
 209             u = u*z + u0[i];
 210             v = v*z + v0[i];



 211         }
 212         return(u/v + tpi*(j0l(x)*logl(x)));

 213 }
 214 
 215 static const GENERIC pr0[12] = {        /* [16 -- inf] */
 216    9.999999999999999999999999999999999997515e-0001L,
 217    1.065981615377273376425365823967550598358e+0003L,
 218    4.390991200927588978306374718984240719130e+0005L,
 219    9.072086218607986711847069407339321363103e+0007L,
 220    1.022552886177375367408408501046461671528e+0010L,
 221    6.420766912243658241570635854089597269031e+0011L,
 222    2.206451725126933913591080211081242266908e+0013L,
 223    3.928369596816895077363705478743346298368e+0014L,
 224    3.258159928874124597286701119721482876596e+0015L,
 225    1.025715808134188978860679130140685101348e+0016L,
 226    7.537170874795721255796001687024031280685e+0015L,
 227   -1.579413901450157332307745586004207687796e+0014L,
 228 };

 229 static const GENERIC ps0[11] = {
 230    1.0e0L,
 231    1.066051927877273376425365823967550512687e+0003L,
 232    4.391739647168381592399173804329266353038e+0005L,
 233    9.075162261801343671805658294123888867884e+0007L,
 234    1.023186118519904751819581912075985995058e+0010L,
 235    6.427861860414223746340515376512730275061e+0011L,
 236    2.210861503237823589735481303627993406235e+0013L,
 237    3.943247335784292905915956840901818177989e+0014L,
 238    3.283720976777545142150200110647270004481e+0015L,
 239    1.045346918812754048903645641538728986759e+0016L,
 240    8.043455468065618900750599584291193680463e+0015L,
 241 };

 242 static const GENERIC pr1[12] = {        /* [8 -- 16] */
 243    9.999999999999999999999784422701108683618e-0001L,
 244    6.796098532948334207755488692777907062894e+0002L,
 245    1.840036112605722168824530758797169836042e+0005L,
 246    2.598490483191916637264894340635847598122e+0007L,
 247    2.105774863242707025525730249472054578523e+0009L,
 248    1.015822044230542426666314997796944979959e+0011L,
 249    2.931557457008110436764077699944189071875e+0012L,
 250    4.962885121125457633655259224179322808824e+0013L,
 251    4.705424055148223269155430598563351566279e+0014L,
 252    2.294439854910747229152056080910427001110e+0015L,
 253    4.905531843137486691500950019322475458629e+0015L,
 254    3.187543169710339218793442542845735994565e+0015L,
 255 };

 256 static const GENERIC ps1[14] = {
 257    1.0e0L,
 258    6.796801657948334207754571576066758180288e+0002L,
 259    1.840512891201300567325421059826676366447e+0005L,
 260    2.599777028312918975306252167127695075221e+0007L,
 261    2.107582572771047636846811284634244892537e+0009L,
 262    1.017275794694156108975782763889979940348e+0011L,
 263    2.938487645192463845428059755454762316011e+0012L,
 264    4.982512164735557054521042916182317924466e+0013L,
 265    4.737639900153703274792677468264564361437e+0014L,
 266    2.323398719123742743524249528275097100646e+0015L,
 267    5.033419107069210577868909797896984419391e+0015L,
 268    3.409036105931068609601317076759804716059e+0015L,
 269    7.505655364352679737585745147753521662166e+0013L,
 270   -9.976837153983688250780198248297109118313e+0012L,
 271 };

 272 static const GENERIC pr2[12] = {        /* [5 -- 8 ] */
 273    9.999999999999999937857236789277366320220e-0001L,
 274    3.692848765268649571651602420376358849214e+0002L,
 275    5.373022067535476576926715900057760985410e+0004L,
 276    4.038738891191314969971504035057219430725e+0006L,
 277    1.728285706306940523397385566659762646999e+0008L,
 278    4.375400819645889911158688737206054788534e+0009L,
 279    6.598950418204912408375591217782088567076e+0010L,
 280    5.827182039183238492480275401520072793783e+0011L,
 281    2.884222642913492390887572414999490975844e+0012L,
 282    7.373278873797767721932837830628688632775e+0012L,
 283    8.338295457568973761205077964397969230489e+0012L,
 284    2.911383183467288345772308817209806922143e+0012L,
 285 };

 286 static const GENERIC ps2[14] = {
 287    1.0e0L,
 288    3.693551890268649477288896267171993213102e+0002L,
 289    5.375607880998361502474715133828068514297e+0004L,
 290    4.042477764024108249744998862572786367328e+0006L,
 291    1.731069838737016956685839588670132939513e+0008L,
 292    4.387147674049898778738226585935491417728e+0009L,
 293    6.628058659620653765349556940567715258165e+0010L,
 294    5.869659904164177740471685856367322160664e+0011L,
 295    2.919839445622817017058977559638969436383e+0012L,
 296    7.535314897696671402628203718612309253907e+0012L,
 297    8.696355561452933775773309859748610658935e+0012L,
 298    3.216155103141537221173601557697083216257e+0012L,
 299    4.756857081068942248246880159213789086363e+0010L,
 300   -3.496356619666608032231074866481472824067e+0009L,
 301 };

 302 static const GENERIC pr3[13] = {        /* [3.5 -- 5 ] */
 303    9.999999999999916693107285612398196588247e-0001L,
 304    2.263975921282917721194425320484974336945e+0002L,
 305    1.994358386744245848889492762781484199966e+0004L,
 306    8.980067458430542243559962493831661323168e+0005L,
 307    2.282213787521372663705567756420087553508e+0007L,
 308    3.409784374889063618250288699908375135923e+0008L,
 309    3.024380857401448589254343517589811711108e+0009L,
 310    1.571110368046740246895071721443082286379e+0010L,
 311    4.603187020243604632153685300463160593768e+0010L,
 312    7.087196453409712719449549280664058793403e+0010L,
 313    5.046196021776346356803687409644239065041e+0010L,
 314    1.287758439080165765709154276618854799932e+0010L,
 315    5.900679773415023433787846658096813590784e+0008L,
 316 };

 317 static const GENERIC ps3[13] = {
 318    1.0e0L,
 319    2.264679046282855061328604619231774747116e+0002L,
 320    1.995939523988944553755653255389812103448e+0004L,
 321    8.993853144706348727038389967490183236820e+0005L,
 322    2.288326099634588843906989983704795468773e+0007L,
 323    3.424967100255240885169240956804790118282e+0008L,
 324    3.046311797972463991368023759640028910016e+0009L,
 325    1.589614961932826812790222479700797224003e+0010L,
 326    4.692406624527744816497089139325073939927e+0010L,
 327    7.320486495902008912866462849073108323948e+0010L,
 328    5.345945972828978289935309597742981360994e+0010L,
 329    1.444033091910423754121309915092247171008e+0010L,
 330    7.987714685115314668378957273824383610525e+0008L,
 331 };

 332 static const GENERIC pr4[13] = {        /* [2.5, 3.5] */
 333    9.999999999986736677961118722747757712260e-0001L,
 334    1.453824980703800559037873123568378845663e+0002L,
 335    8.097327216430682288267610447006508661032e+0003L,
 336    2.273847252038264370231169686380192662135e+0005L,
 337    3.561056728046211111354759998976985449622e+0006L,
 338    3.244933588800096378434627029369680378599e+0007L,
 339    1.740112392860717950376210038908476792588e+0008L,
 340    5.426170187455893285197878563881579269524e+0008L,
 341    9.490107486454362321004377336020526281371e+0008L,
 342    8.688872439428470049801714121070005313806e+0008L,
 343    3.673315853166437222811910656900123215515e+0008L,
 344    5.577770470359303305164877446339693270239e+0007L,
 345    1.540438642031689641308197880181291865714e+0006L,
 346 };

 347 static const GENERIC ps4[13] = {        /* [2.5, 3.5] */
 348    1.0e0L,
 349    1.454528105698159439773035951959131799816e+0002L,
 350    8.107442215200392397172179900434987859618e+0003L,
 351    2.279390393778242887574177096606328994140e+0005L,
 352    3.576251625592252008424781111770934135844e+0006L,
 353    3.267909499056932631405942058670933813863e+0007L,
 354    1.760021515330805537499778238099704648805e+0008L,
 355    5.525553787667353981242060222587465726729e+0008L,
 356    9.769870295912820457889384082671269328511e+0008L,
 357    9.110582071004774279226905629624018008454e+0008L,
 358    3.981857678621955599371967680343918454345e+0008L,
 359    6.482404686230769399073192961667697036706e+0007L,
 360    2.210046943095878402443535460329391782298e+0006L,
 361 };

 362 static const GENERIC pr5[13] = {        /* [1.777..., 2.5] */
 363    9.999999999114986107951817871144655880699e-0001L,
 364    9.252583736048588342568344570315435947614e+0001L,
 365    3.218726757856078715214631502407386264637e+0003L,
 366    5.554009964621111656479588505862577040831e+0004L,
 367    5.269993115643664338253196944523510290175e+0005L,
 368    2.874613773778430691192912190618220544575e+0006L,
 369    9.133538151103658353874146919613442436035e+0006L,
 370    1.673067041410338922825193013077354249193e+0007L,
 371    1.706913873848398011744790289200151840498e+0007L,
 372    9.067766583853288534551600235576747618679e+0006L,
 373    2.216746733457884568532695355036338655872e+0006L,
 374    1.945753880802872541235703812722344514405e+0005L,
 375    3.132374412921948071539195638885330951749e+0003L,
 376 };

 377 static const GENERIC ps5[13] = {        /* [1.777..., 2.5] */
 378    1.0e0L,
 379    9.259614983862181118883831670990340052982e+0001L,
 380    3.225125275462903384842124075132609290304e+0003L,
 381    5.575705362829101545292760055941855246492e+0004L,
 382    5.306049863037087855496170121958448492522e+0005L,
 383    2.907060758873509564309729903109018597215e+0006L,
 384    9.298059206584995898298257827131208539289e+0006L,
 385    1.720391071006963176836108026556547062980e+0007L,
 386    1.782614812922865190479394509487941920612e+0007L,
 387    9.708016389605273153536452032839879950155e+0006L,
 388    2.476495084688170096480215640962175140027e+0006L,
 389    2.363200660365585759668077790194604917187e+0005L,
 390    4.803239569848196077121203575704356936731e+0003L,
 391 };

 392 static const GENERIC pr6[13] = {        /* [1.28, 1.777...] */
 393    9.999999969777095495998606925524322559556e-0001L,
 394    5.825486719466194430503283824096872219216e+0001L,
 395    1.248155491637757281915184824965379905380e+0003L,
 396    1.302093199842358609321338417071710477615e+0004L,
 397    7.353835804186292782840961999810543016039e+0004L,
 398    2.356471661113686180549195092555751341757e+0005L,
 399    4.350553267429009581632987060942780847101e+0005L,
 400    4.588762661876600638719159826652389418235e+0005L,
 401    2.675796398548523436544221045225290128611e+0005L,
 402    8.077649557108971388298292919988449940464e+0004L,
 403    1.117640459221306873519068741664054573776e+0004L,
 404    5.544400072396814695175787511557757885585e+0002L,
 405    5.072550541191480498431289089905822910718e+0000L,
 406 };

 407 static const GENERIC ps6[13] = {        /* [1.28, 1.777...] */
 408    1.0e0L,
 409    5.832517925357165050639075848183613063291e+0001L,
 410    1.252144364743592128171256104364976466898e+0003L,
 411    1.310300234342216813579118022415585740772e+0004L,
 412    7.434667697093812197817292154032863632923e+0004L,
 413    2.398706595587719165726469002404004614711e+0005L,
 414    4.472737517625103157004869372427480602511e+0005L,
 415    4.786313523337761975294171429067037723611e+0005L,
 416    2.851161872872731228472536061865365370192e+0005L,
 417    8.891648269899148412331918021801385815586e+0004L,
 418    1.297097489535351517572978123584751042287e+0004L,
 419    7.096761640545975756202184143400469812618e+0002L,
 420    8.378049338590233325977702401733340820351e+0000L,
 421 };

 422 static const GENERIC sixteen = 16.0L;
 423 static const GENERIC huge    = 1.0e30L;
 424 
 425 static GENERIC pzero(x)
 426 GENERIC x;
 427 {
 428         GENERIC s,r,t,z;
 429         int i;
 430         if (x>huge) return one;
 431         t = one/x; z = t*t;
 432         if (x>sixteen) {
 433             r = z*pr0[11]+pr0[10]; s = ps0[10];
 434             for (i=9;i>=0;i--) {
 435                 r = z*r + pr0[i];
 436                 s = z*s + ps0[i];







 437             }
 438         } else if (x > eight) {
 439             r = pr1[11]; s = ps1[11]+z*(ps1[12]+z*ps1[13]);
 440             for (i=10;i>=0;i--) {
 441                 r = z*r + pr1[i];
 442                 s = z*s + ps1[i];


 443             }
 444         } else if (x > five) {       /* x > 5.0 */
 445             r = pr2[11]; s = ps2[11]+z*(ps2[12]+z*ps2[13]);
 446             for (i=10;i>=0;i--) {
 447                 r = z*r + pr2[i];
 448                 s = z*s + ps2[i];
 449             }
 450         } else if (x>3.5L) {
 451             r = pr3[12]; s = ps3[12];
 452             for (i=11;i>=0;i--) {
 453                 r = z*r + pr3[i];
 454                 s = z*s + ps3[i];
 455             }
 456         } else if (x>2.5L) {
 457             r = pr4[12]; s = ps4[12];
 458             for (i=11;i>=0;i--) {
 459                 r = z*r + pr4[i];
 460                 s = z*s + ps4[i];
 461             }
 462         } else if (x> (1.0L/0.5625L)) {
 463             r = pr5[12]; s = ps5[12];
 464             for (i=11;i>=0;i--) {
 465                 r = z*r + pr5[i];
 466                 s = z*s + ps5[i];








 467             }
 468         } else {        /* assume x > 1.28 */
 469             r = pr6[12]; s = ps6[12];
 470             for (i=11;i>=0;i--) {
 471                 r = z*r + pr6[i];
 472                 s = z*s + ps6[i];


 473             }
 474         }
 475         return r/s;
 476 }
 477 
 478 


 479 static const GENERIC qr0[12] = {        /* [16, inf] */
 480   -1.249999999999999999999999999999999972972e-0001L,
 481   -1.425179595545670577414395762503991596897e+0002L,
 482   -6.312499645625970845534460257936222407219e+0004L,
 483   -1.411374326457208384315121243698814446848e+0007L,
 484   -1.735034212758873581410984757860787252842e+0009L,
 485   -1.199777647512789489421826342485055280680e+0011L,
 486   -4.596025334081655714499860409699100373644e+0012L,
 487   -9.262525628201284107792924477031653399187e+0013L,
 488   -8.858394728685039245344398842180662867639e+0014L,
 489   -3.267527953687534887623740622709505972113e+0015L,
 490   -2.664222971186311967587129347029450062019e+0015L,
 491    3.442464060723987869585180095344504100204e+0014L,
 492 };

 493 static const GENERIC qs0[11] = {
 494    1.0e0L,
 495    1.140729613936536461931516610003185687881e+0003L,
 496    5.056665510442299351009198186490085803580e+0005L,
 497    1.132041763825642787943941650522718199115e+0008L,
 498    1.394570111872581606392620678214246479767e+0010L,
 499    9.677945218152264789534431079563744378421e+0011L,
 500    3.731140327851536828225143058896348502096e+0013L,
 501    7.612785951064869291722846681020881676410e+0014L,
 502    7.476077016406764891730191004811863975940e+0015L,
 503    2.951246482613592035421503427100393831709e+0016L,
 504    3.108361803691811711136854587074302034901e+0016L,
 505 };

 506 static const GENERIC qr1[12] = {        /* [8, 16 ] */
 507   -1.249999999999999999997949010383433818157e-0001L,
 508   -9.051215166393822640636752244895124126934e+0001L,
 509   -2.620782703428148837671179031904208303947e+0004L,
 510   -3.975571261553504457766177974508785790884e+0006L,
 511   -3.479029330759311306270072218074074994090e+0008L,
 512   -1.823955008124268573036216746186239829089e+0010L,
 513   -5.765932697111801375765156029221568664435e+0011L,
 514   -1.079843680798742592954002192417934779114e+0013L,
 515   -1.146893630504592739082205764611581332897e+0014L,
 516   -6.367016059683898464936104447282880704182e+0014L,
 517   -1.583109041961213490464459111903484209098e+0015L,
 518   -1.230149555764242473103128650135795639412e+0015L,
 519 };

 520 static const GENERIC qs1[14] = {
 521    1.0e0L,
 522    7.246831508115058112438579847778014458432e+0002L,
 523    2.100854184439168518399383786306927037611e+0005L,
 524    3.192636418837951507430188285940994235122e+0007L,
 525    2.801558443383354674538443461124434216152e+0009L,
 526    1.475026997664373739293483927250653467487e+0011L,
 527    4.694486824913954608552363821799927145318e+0012L,
 528    8.890350100919200250838438709601547334021e+0013L,
 529    9.626844429082905144874701068760469752067e+0014L,
 530    5.541110744600460773528263862687521642140e+0015L,
 531    1.486500494789452556727470329232123096563e+0016L,
 532    1.415840104845959400365430773732093899210e+0016L,
 533    1.780866095241517418081312567239682336483e+0015L,
 534   -2.359230917384889357887631544079990129494e+0014L,
 535 };

 536 static const GENERIC qr2[12] = {        /* [5, 8] */
 537   -1.249999999999999531937744362527772181614e-0001L,
 538   -4.944373897356969774839375977239241573966e+0001L,
 539   -7.728449175433465285314261650078450473909e+0003L,
 540   -6.262574329612752346336901434651220705903e+0005L,
 541   -2.900948220220943306027235217424380672732e+0007L,
 542   -7.988719647634192770463917157562874119535e+0008L,
 543   -1.318228171927181389547760026626357012375e+0010L,
 544   -1.282439773983029245309263271945424928196e+0011L,
 545   -7.050925570827818040186149940257918845138e+0011L,
 546   -2.021751882573871990004205616874202684429e+0012L,
 547   -2.592939962400668552384333900573812635658e+0012L,
 548   -1.038267109518891262840601514932972850326e+0012L,
 549 };

 550 static const GENERIC qs2[14] = {
 551    1.0e0L,
 552    3.961358492885570003202784022894248952116e+0002L,
 553    6.205788738864701882828752634586510926968e+0004L,
 554    5.045715603932670286550673813011764406749e+0006L,
 555    2.349248611362658323353343389430968751429e+0008L,
 556    6.520244524415828635917683553721880063911e+0009L,
 557    1.089111211223507719337067159886281887722e+0011L,
 558    1.080406000905359867958779409414903018610e+0012L,
 559    6.135645280895514703514154680623769562148e+0012L,
 560    1.862433040246625874245867151368643668215e+0013L,
 561    2.667780805786648888840777888702193708994e+0013L,
 562    1.394401107289087774765300711809313112824e+0013L,
 563    1.093247500616320375562898297156722445484e+0012L,
 564   -7.228875530378928722826604216491493780775e+0010L,
 565 };

 566 static const GENERIC qr3[13] = {        /* [3.5 5] */
 567   -1.249999999999473067748420379578481661075e-0001L,
 568   -3.044549048635289351913574324803250977998e+0001L,
 569   -2.890081140649769078496693003524681440869e+0003L,
 570   -1.404922456817202235879343275330529107684e+0005L,
 571   -3.862746614385573443518177403617349281869e+0006L,
 572   -6.257517309110249049201133708911155047689e+0007L,
 573   -6.031451330920839916987079782727323477520e+0008L,
 574   -3.411542405173830611454025765755854382346e+0009L,
 575   -1.089392478149726672133014498723021526099e+0010L,
 576   -1.824934078420210941290140903415956782726e+0010L,
 577   -1.400780278304358710423481070486939531139e+0010L,
 578   -3.716484136064917363926635716743771092093e+0009L,
 579   -1.397591075296425529970434890954904331580e+0008L,
 580 };

 581 static const GENERIC qs3[13] = {
 582    1.0e0L,
 583    2.441498613904962049391000187014945858042e+0002L,
 584    2.326188882072370711500164222341514337043e+0004L,
 585    1.137138213121231338494977104659239578165e+0006L,
 586    3.152918070735662728722998452605364253517e+0007L,
 587    5.172877993426507259314270488444013595108e+0008L,
 588    5.083086439731669807455961078856470774115e+0009L,
 589    2.961842732066434123119325521139476909941e+0010L,
 590    9.912185866862440735829781856081353151390e+0010L,
 591    1.793560561251622234430564181567297983598e+0011L,
 592    1.577090119341228122525265108497940403073e+0011L,
 593    5.509910306780166194333889999985463681636e+0010L,
 594    4.761691134078874491202320181517936758141e+0009L,
 595 };

 596 static const GENERIC qr4[13] = {        /* [2.5 3.5] */
 597   -1.249999999928567734339745043490705340835e-0001L,
 598   -1.967201748731419063051601624435565528481e+0001L,
 599   -1.186329146714562236407099740615528170707e+0003L,
 600   -3.607736959222941810356301491152457934060e+0004L,
 601   -6.119200717978104904932828468575194267125e+0005L,
 602   -6.037847781158358226670305078652205586384e+0006L,
 603   -3.503558153336140359700536720393565984740e+0007L,
 604   -1.180196478268225718757218523746787309773e+0008L,
 605   -2.221860232085134915841426363505169680528e+0008L,
 606   -2.173372505452747585296176761701746236760e+0008L,
 607   -9.649364865061237558517730539506568013963e+0007L,
 608   -1.465429227847933034546039640094862650385e+0007L,
 609   -3.083003197920262085170581866246663380607e+0005L,
 610 };

 611 static const GENERIC qs4[13] = {        /* [2.5 3.5] */
 612    1.0e0L,
 613    1.579620773732259142752614142139986854055e+0002L,
 614    9.581372220329138733203879503753685054968e+0003L,
 615    2.939598672379108095776114131010825885308e+0005L,
 616    5.052183049314542218630341818692588448168e+0006L,
 617    5.083497695595206639433839326338971980149e+0007L,
 618    3.036385361800553388049719014005099206516e+0008L,
 619    1.067826481452753409910563785161661492137e+0009L,
 620    2.145644125557118044720741775125319669272e+0009L,
 621    2.324115615959719949363946673491552216799e+0009L,
 622    1.223262962112070757966959855619847011146e+0009L,
 623    2.569765553318495423738478585947110270709e+0008L,
 624    1.354744744299227127897905787732636565504e+0007L,
 625 };

 626 static const GENERIC qr5[13] = {        /* [1.777.., 2.5] */
 627   -1.249999995936639697637680428174576069971e-0001L,
 628   -1.260846055371311453485891923426489068315e+0001L,
 629   -4.772398467544467480801174330290141578895e+0002L,
 630   -8.939852599990298486613760833996490599724e+0003L,
 631   -9.184070787149542050979542226446134243197e+0004L,
 632   -5.406038945018274458362637897739280435171e+0005L,
 633   -1.845896544705190261018653728678171084418e+0006L,
 634   -3.613616990680809501878667570653308071547e+0006L,
 635   -3.908782978135693252252557720414348623779e+0006L,
 636   -2.173711022517323927109138170588442768176e+0006L,
 637   -5.431253130679918485836408549007856244495e+0005L,
 638   -4.591098546452684510082591587275940765959e+0004L,
 639   -5.244711364168207806835520057792229646578e+0002L,
 640 };

 641 static const GENERIC qs5[13] = {        /* [1.777.., 2.5] */
 642    1.0e0L,
 643    1.014536210851290878350892750972474861447e+0002L,
 644    3.875547510687135314064434160096139681076e+0003L,
 645    7.361913122670079814955259281995617732580e+0004L,
 646    7.720288944218771126581086539585529314636e+0005L,
 647    4.681529554446752496404431433608306558038e+0006L,
 648    1.667882621940503925455031252308367745820e+0007L,
 649    3.469403153761399881888272620855305156241e+0007L,
 650    4.096992047964210711867089384719947863019e+0007L,
 651    2.596804755829217449311530735959560630554e+0007L,
 652    7.983933774697889238154465064019410763845e+0006L,
 653    9.818133816979900819087242425280757938152e+0005L,
 654    3.061083930868694396013541535670745443560e+0004L,
 655 };
 656 
 657 static const GENERIC qr6[13] = {        /* [1.28, 1.777..] */
 658   -1.249999881577289001807137282824929082771e-0001L,
 659   -7.998273510053110759610810594119533619282e+0000L,
 660   -1.872481955335172543369089617771565632719e+0002L,
 661   -2.122116786726300805079874003303799646812e+0003L,
 662   -1.293850285839529282503178263484773478457e+0004L,
 663   -4.445024742266316181033354192262529356093e+0004L,
 664   -8.730161378334357767668344467356505347070e+0004L,
 665   -9.706222895172078442801444972505315054736e+0004L,
 666   -5.896325518259858270165531513618195321041e+0004L,
 667   -1.823172034368108822276420827074668832233e+0004L,
 668   -2.509304178635055926638833040337472387175e+0003L,
 669   -1.156608965715779237316769828941729964099e+0002L,
 670   -7.028005789650731396887346826397785210442e-0001L,
 671 };

 672 static const GENERIC qs6[13] = {        /* [1.28, 1.777..] */
 673    1.0e0L,
 674    6.457211085058064845601261321277721075900e+0001L,
 675    1.534005216588011210342824555136008682950e+0003L,
 676    1.777217999176441782593357660462379097171e+0004L,
 677    1.118372652642469468091084810263231199696e+0005L,
 678    4.015242433858461813142365748386473605294e+0005L,
 679    8.377081045517098645448616514388280497673e+0005L,
 680    1.011495020008010352575398009604164287337e+0006L,
 681    6.886722075290430568652227875200208955970e+0005L,
 682    2.504735189948021472047157148613171956537e+0005L,
 683    4.408138920171044846941001844352009817062e+0004L,
 684    3.105572178072115145673058722853640854884e+0003L,
 685    5.588294821118916113437396504573817033678e+0001L,
 686 };
 687 static GENERIC qzero(x)
 688 GENERIC x;

 689 {
 690         GENERIC s,r,t,z;
 691         int i;
 692         if (x>huge) return -0.125L/x;
 693         t = one/x; z = t*t;
 694         if (x>sixteen) {
 695             r = z*qr0[11]+qr0[10]; s = qs0[10];
 696             for (i=9;i>=0;i--) {
 697                 r = z*r + qr0[i];
 698                 s = z*s + qs0[i];
 699             }
 700         } else if (x>eight) {
 701             r = qr1[11]; s = qs1[11]+z*(qs1[12]+z*qs1[13]);
 702             for (i=10;i>=0;i--) {
 703                 r = z*r + qr1[i];
 704                 s = z*s + qs1[i];
 705             }
 706         } else if (x>five) {  /* assume x > 5.0 */
 707             r = qr2[11]; s = qs2[11]+z*(qs2[12]+z*qs2[13]);
 708             for (i=10;i>=0;i--) {
 709                 r = z*r + qr2[i];
 710                 s = z*s + qs2[i];
 711             }
 712         } else if (x>3.5L) {
 713             r = qr3[12]; s = qs3[12];
 714             for (i=11;i>=0;i--) {
 715                 r = z*r + qr3[i];
 716                 s = z*s + qs3[i];
 717             }
 718         } else if (x>2.5L) {
 719             r = qr4[12]; s = qs4[12];
 720             for (i=11;i>=0;i--) {
 721                 r = z*r + qr4[i];
 722                 s = z*s + qs4[i];
 723             }
 724         } else if (x> (1.0L/0.5625L)) {
 725             r = qr5[12]; s = qs5[12];
 726             for (i=11;i>=0;i--) {
 727                 r = z*r + qr5[i];
 728                 s = z*s + qs5[i];

















 729             }
 730         } else {        /* assume x > 1.28 */
 731             r = qr6[12]; s = qs6[12];
 732             for (i=11;i>=0;i--) {
 733                 r = z*r + qr6[i];
 734                 s = z*s + qs6[i];


 735             }
 736         }
 737         return t*(r/s);

 738 }


   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21 
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 
  26 /*
  27  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  28  * Use is subject to license terms.
  29  */
  30 
  31 /*
  32  * Floating point Bessel's function of the first and second kinds
  33  * of order zero: j0(x),y0(x);
  34  *
  35  * Special cases:
  36  *      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
  37  *      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
  38  */
  39 
  40 #pragma weak __j0l = j0l
  41 #pragma weak __y0l = y0l
  42 
  43 #include "libm.h"
  44 #include "longdouble.h"
  45 
  46 #define GENERIC long double
  47 
  48 static const GENERIC zero = 0.0L,
  49         small = 1.0e-9L,
  50         tiny = 1.0e-38L,
  51         one = 1.0L,
  52         five = 5.0L,
  53         eight = 8.0L,
  54         invsqrtpi = 5.641895835477562869480794515607725858441e-0001L,
  55         tpi = 0.636619772367581343075535053490057448L;
  56 
  57 static GENERIC pzero(GENERIC);
  58 static GENERIC qzero(GENERIC);
  59 
  60 static GENERIC r0[7] = {
  61         -2.499999999999999999999999999999998934492e-0001L,
  62         1.272657927360049786327618451133763714880e-0002L,
  63         -2.694499763712963276900636693400659600898e-0004L,
  64         2.724877475058977576903234070919616447883e-0006L,
  65         -1.432617103214330236967477495393076320281e-0008L,
  66         3.823248804080079168706683540513792224471e-0011L,
  67         -4.183174277567983647337568504286313665065e-0014L,
  68 };
  69 
  70 static GENERIC s0[7] = {
  71         1.0e0L,
  72         1.159368290559800854689526195462884666395e-0002L,
  73         6.629397597394973383009743876169946772559e-0005L,
  74         2.426779981394054406305431142501735094340e-0007L,
  75         6.097663491248511069094400469635449749883e-0010L,
  76         1.017019133340929220238747413216052224036e-0012L,
  77         9.012593179306197579518374581969371278481e-0016L,
  78 };
  79 
  80 GENERIC
  81 j0l(GENERIC x)
  82 {
  83         GENERIC z, s, c, ss, cc, r, u, v;
  84         int i;
  85 
  86         if (isnanl(x))
  87                 return (x + x);
  88 
  89         x = fabsl(x);
  90 
  91         if (x > 1.28L) {
  92                 if (!finitel(x))
  93                         return (zero);
  94 
  95                 s = sinl(x);
  96                 c = cosl(x);
  97 
  98                 /* BEGIN CSTYLED */
  99                 /*
 100                  * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
 101                  * where x0 = x-pi/4
 102                  *      Better formula:
 103                  *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
 104                  *                      = 1/sqrt(2) * (cos(x) + sin(x))
 105                  *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
 106                  *                      = 1/sqrt(2) * (sin(x) - cos(x))
 107                  * To avoid cancellation, use
 108                  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 109                  * to compute the worse one.
 110                  */
 111                 /* END CSTYLED */
 112                 if (x > 1.0e2450L) { /* x+x may overflow */
 113                         ss = s - c;
 114                         cc = s + c;
 115                 } else if (signbitl(s) != signbitl(c)) {
 116                         ss = s - c;
 117                         cc = -cosl(x + x) / ss;
 118                 } else {
 119                         cc = s + c;
 120                         ss = -cosl(x + x) / cc;
 121                 }
 122 
 123                 /*
 124                  * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 125                  * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 126                  */
 127                 if (x > 1.0e120L)
 128                         return ((invsqrtpi * cc) / sqrtl(x));
 129 
 130                 u = pzero(x);
 131                 v = qzero(x);
 132                 return (invsqrtpi * (u * cc - v * ss) / sqrtl(x));
 133         }
 134 
 135         if (x <= small) {
 136                 if (x <= tiny)
 137                         return (one - x);
 138                 else
 139                         return (one - x * x * 0.25L);
 140         }
 141 
 142         z = x * x;
 143         r = r0[6];
 144         s = s0[6];
 145 
 146         for (i = 5; i >= 0; i--) {
 147                 r = r * z + r0[i];
 148                 s = s * z + s0[i];
 149         }
 150 
 151         return (one + z * (r / s));
 152 }
 153 
 154 static const GENERIC u0[8] = {
 155         -7.380429510868722527434392794848301631220e-0002L,
 156         1.766855559625940791857536949301981816513e-0001L,
 157         -1.386470722701047923235553251240162839408e-0002L,
 158         3.520149242724811578636970811631224862615e-0004L,
 159         -3.978599663243790049853642275624951870025e-0006L,
 160         2.228801153263957224547222556806915479763e-0008L,
 161         -6.121246764298785018658597179498837316177e-0011L,
 162         6.677103629722678833475965810525587396596e-0014L,
 163 };
 164 
 165 static const GENERIC v0[8] = {
 166         1.0e0L,
 167         1.247164416539111311571676766127767127970e-0002L,
 168         7.829144749639791500052900281489367443576e-0005L,
 169         3.247126540422245330511218321013360336606e-0007L,
 170         9.750516724789499678567062572549568447869e-0010L,
 171         2.156713223173591212250543390258458098776e-0012L,
 172         3.322169561597890004231482431236452752624e-0015L,
 173         2.821213295314000924252226486305726805093e-0018L,
 174 };
 175 
 176 GENERIC
 177 y0l(GENERIC x)
 178 {
 179         GENERIC z, s, c, ss, cc, u, v;
 180         int i;
 181         volatile GENERIC d;
 182 
 183         if (isnanl(x))
 184                 return (x + x);
 185 
 186         if (x <= zero) {
 187                 if (x == zero)
 188                         d = -one / (x - x);
 189                 else
 190                         d = zero / (x - x);
 191         }
 192 
 193 #ifdef lint
 194         d = d;
 195 #endif
 196 
 197         if (x > 1.28L) {
 198                 if (!finitel(x))
 199                         return (zero);
 200 
 201                 s = sinl(x);
 202                 c = cosl(x);
 203 
 204                 /* BEGIN CSTYLED */
 205                 /*
 206                  * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
 207                  * where x0 = x-pi/4
 208                  *      Better formula:
 209                  *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
 210                  *                      = 1/sqrt(2) * (cos(x) + sin(x))
 211                  *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
 212                  *                      = 1/sqrt(2) * (sin(x) - cos(x))
 213                  * To avoid cancellation, use
 214                  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 215                  * to compute the worse one.
 216                  */
 217                 /* END CSTYLED */
 218                 if (x > 1.0e2450L) { /* x+x may overflow */
 219                         ss = s - c;
 220                         cc = s + c;
 221                 } else if (signbitl(s) != signbitl(c)) {
 222                         ss = s - c;
 223                         cc = -cosl(x + x) / ss;
 224                 } else {
 225                         cc = s + c;
 226                         ss = -cosl(x + x) / cc;
 227                 }
 228 
 229                 /*
 230                  * j0(x) = 1/sqrt(pi*x) * (P(0,x)*cc - Q(0,x)*ss)
 231                  * y0(x) = 1/sqrt(pi*x) * (P(0,x)*ss + Q(0,x)*cc)
 232                  */
 233                 if (x > 1.0e120L)
 234                         return ((invsqrtpi * ss) / sqrtl(x));
 235 
 236                 return (invsqrtpi * (pzero(x) * ss + qzero(x) * cc) / sqrtl(x));
 237         }
 238 
 239         if (x <= tiny)
 240                 return (u0[0] + tpi * logl(x));
 241 
 242         z = x * x;
 243         u = u0[7];
 244         v = v0[7];
 245 
 246         for (i = 6; i >= 0; i--) {
 247                 u = u * z + u0[i];
 248                 v = v * z + v0[i];
 249         }
 250 
 251         return (u / v + tpi * (j0l(x) * logl(x)));
 252 }
 253 
 254 static const GENERIC pr0[12] = {        /* [16 -- inf] */
 255         9.999999999999999999999999999999999997515e-0001L,
 256         1.065981615377273376425365823967550598358e+0003L,
 257         4.390991200927588978306374718984240719130e+0005L,
 258         9.072086218607986711847069407339321363103e+0007L,
 259         1.022552886177375367408408501046461671528e+0010L,
 260         6.420766912243658241570635854089597269031e+0011L,
 261         2.206451725126933913591080211081242266908e+0013L,
 262         3.928369596816895077363705478743346298368e+0014L,
 263         3.258159928874124597286701119721482876596e+0015L,
 264         1.025715808134188978860679130140685101348e+0016L,
 265         7.537170874795721255796001687024031280685e+0015L,
 266         -1.579413901450157332307745586004207687796e+0014L,
 267 };
 268 
 269 static const GENERIC ps0[11] = {
 270         1.0e0L,
 271         1.066051927877273376425365823967550512687e+0003L,
 272         4.391739647168381592399173804329266353038e+0005L,
 273         9.075162261801343671805658294123888867884e+0007L,
 274         1.023186118519904751819581912075985995058e+0010L,
 275         6.427861860414223746340515376512730275061e+0011L,
 276         2.210861503237823589735481303627993406235e+0013L,
 277         3.943247335784292905915956840901818177989e+0014L,
 278         3.283720976777545142150200110647270004481e+0015L,
 279         1.045346918812754048903645641538728986759e+0016L,
 280         8.043455468065618900750599584291193680463e+0015L,
 281 };
 282 
 283 static const GENERIC pr1[12] = {        /* [8 -- 16] */
 284         9.999999999999999999999784422701108683618e-0001L,
 285         6.796098532948334207755488692777907062894e+0002L,
 286         1.840036112605722168824530758797169836042e+0005L,
 287         2.598490483191916637264894340635847598122e+0007L,
 288         2.105774863242707025525730249472054578523e+0009L,
 289         1.015822044230542426666314997796944979959e+0011L,
 290         2.931557457008110436764077699944189071875e+0012L,
 291         4.962885121125457633655259224179322808824e+0013L,
 292         4.705424055148223269155430598563351566279e+0014L,
 293         2.294439854910747229152056080910427001110e+0015L,
 294         4.905531843137486691500950019322475458629e+0015L,
 295         3.187543169710339218793442542845735994565e+0015L,
 296 };
 297 
 298 static const GENERIC ps1[14] = {
 299         1.0e0L,
 300         6.796801657948334207754571576066758180288e+0002L,
 301         1.840512891201300567325421059826676366447e+0005L,
 302         2.599777028312918975306252167127695075221e+0007L,
 303         2.107582572771047636846811284634244892537e+0009L,
 304         1.017275794694156108975782763889979940348e+0011L,
 305         2.938487645192463845428059755454762316011e+0012L,
 306         4.982512164735557054521042916182317924466e+0013L,
 307         4.737639900153703274792677468264564361437e+0014L,
 308         2.323398719123742743524249528275097100646e+0015L,
 309         5.033419107069210577868909797896984419391e+0015L,
 310         3.409036105931068609601317076759804716059e+0015L,
 311         7.505655364352679737585745147753521662166e+0013L,
 312         -9.976837153983688250780198248297109118313e+0012L,
 313 };
 314 
 315 static const GENERIC pr2[12] = {        /* [5 -- 8 ] */
 316         9.999999999999999937857236789277366320220e-0001L,
 317         3.692848765268649571651602420376358849214e+0002L,
 318         5.373022067535476576926715900057760985410e+0004L,
 319         4.038738891191314969971504035057219430725e+0006L,
 320         1.728285706306940523397385566659762646999e+0008L,
 321         4.375400819645889911158688737206054788534e+0009L,
 322         6.598950418204912408375591217782088567076e+0010L,
 323         5.827182039183238492480275401520072793783e+0011L,
 324         2.884222642913492390887572414999490975844e+0012L,
 325         7.373278873797767721932837830628688632775e+0012L,
 326         8.338295457568973761205077964397969230489e+0012L,
 327         2.911383183467288345772308817209806922143e+0012L,
 328 };
 329 
 330 static const GENERIC ps2[14] = {
 331         1.0e0L,
 332         3.693551890268649477288896267171993213102e+0002L,
 333         5.375607880998361502474715133828068514297e+0004L,
 334         4.042477764024108249744998862572786367328e+0006L,
 335         1.731069838737016956685839588670132939513e+0008L,
 336         4.387147674049898778738226585935491417728e+0009L,
 337         6.628058659620653765349556940567715258165e+0010L,
 338         5.869659904164177740471685856367322160664e+0011L,
 339         2.919839445622817017058977559638969436383e+0012L,
 340         7.535314897696671402628203718612309253907e+0012L,
 341         8.696355561452933775773309859748610658935e+0012L,
 342         3.216155103141537221173601557697083216257e+0012L,
 343         4.756857081068942248246880159213789086363e+0010L,
 344         -3.496356619666608032231074866481472824067e+0009L,
 345 };
 346 
 347 static const GENERIC pr3[13] = {        /* [3.5 -- 5 ] */
 348         9.999999999999916693107285612398196588247e-0001L,
 349         2.263975921282917721194425320484974336945e+0002L,
 350         1.994358386744245848889492762781484199966e+0004L,
 351         8.980067458430542243559962493831661323168e+0005L,
 352         2.282213787521372663705567756420087553508e+0007L,
 353         3.409784374889063618250288699908375135923e+0008L,
 354         3.024380857401448589254343517589811711108e+0009L,
 355         1.571110368046740246895071721443082286379e+0010L,
 356         4.603187020243604632153685300463160593768e+0010L,
 357         7.087196453409712719449549280664058793403e+0010L,
 358         5.046196021776346356803687409644239065041e+0010L,
 359         1.287758439080165765709154276618854799932e+0010L,
 360         5.900679773415023433787846658096813590784e+0008L,
 361 };
 362 
 363 static const GENERIC ps3[13] = {
 364         1.0e0L,
 365         2.264679046282855061328604619231774747116e+0002L,
 366         1.995939523988944553755653255389812103448e+0004L,
 367         8.993853144706348727038389967490183236820e+0005L,
 368         2.288326099634588843906989983704795468773e+0007L,
 369         3.424967100255240885169240956804790118282e+0008L,
 370         3.046311797972463991368023759640028910016e+0009L,
 371         1.589614961932826812790222479700797224003e+0010L,
 372         4.692406624527744816497089139325073939927e+0010L,
 373         7.320486495902008912866462849073108323948e+0010L,
 374         5.345945972828978289935309597742981360994e+0010L,
 375         1.444033091910423754121309915092247171008e+0010L,
 376         7.987714685115314668378957273824383610525e+0008L,
 377 };
 378 
 379 static const GENERIC pr4[13] = {        /* [2.5, 3.5] */
 380         9.999999999986736677961118722747757712260e-0001L,
 381         1.453824980703800559037873123568378845663e+0002L,
 382         8.097327216430682288267610447006508661032e+0003L,
 383         2.273847252038264370231169686380192662135e+0005L,
 384         3.561056728046211111354759998976985449622e+0006L,
 385         3.244933588800096378434627029369680378599e+0007L,
 386         1.740112392860717950376210038908476792588e+0008L,
 387         5.426170187455893285197878563881579269524e+0008L,
 388         9.490107486454362321004377336020526281371e+0008L,
 389         8.688872439428470049801714121070005313806e+0008L,
 390         3.673315853166437222811910656900123215515e+0008L,
 391         5.577770470359303305164877446339693270239e+0007L,
 392         1.540438642031689641308197880181291865714e+0006L,
 393 };
 394 
 395 static const GENERIC ps4[13] = {        /* [2.5, 3.5] */
 396         1.0e0L,
 397         1.454528105698159439773035951959131799816e+0002L,
 398         8.107442215200392397172179900434987859618e+0003L,
 399         2.279390393778242887574177096606328994140e+0005L,
 400         3.576251625592252008424781111770934135844e+0006L,
 401         3.267909499056932631405942058670933813863e+0007L,
 402         1.760021515330805537499778238099704648805e+0008L,
 403         5.525553787667353981242060222587465726729e+0008L,
 404         9.769870295912820457889384082671269328511e+0008L,
 405         9.110582071004774279226905629624018008454e+0008L,
 406         3.981857678621955599371967680343918454345e+0008L,
 407         6.482404686230769399073192961667697036706e+0007L,
 408         2.210046943095878402443535460329391782298e+0006L,
 409 };
 410 
 411 static const GENERIC pr5[13] = {        /* [1.777..., 2.5] */
 412         9.999999999114986107951817871144655880699e-0001L,
 413         9.252583736048588342568344570315435947614e+0001L,
 414         3.218726757856078715214631502407386264637e+0003L,
 415         5.554009964621111656479588505862577040831e+0004L,
 416         5.269993115643664338253196944523510290175e+0005L,
 417         2.874613773778430691192912190618220544575e+0006L,
 418         9.133538151103658353874146919613442436035e+0006L,
 419         1.673067041410338922825193013077354249193e+0007L,
 420         1.706913873848398011744790289200151840498e+0007L,
 421         9.067766583853288534551600235576747618679e+0006L,
 422         2.216746733457884568532695355036338655872e+0006L,
 423         1.945753880802872541235703812722344514405e+0005L,
 424         3.132374412921948071539195638885330951749e+0003L,
 425 };
 426 
 427 static const GENERIC ps5[13] = {        /* [1.777..., 2.5] */
 428         1.0e0L,
 429         9.259614983862181118883831670990340052982e+0001L,
 430         3.225125275462903384842124075132609290304e+0003L,
 431         5.575705362829101545292760055941855246492e+0004L,
 432         5.306049863037087855496170121958448492522e+0005L,
 433         2.907060758873509564309729903109018597215e+0006L,
 434         9.298059206584995898298257827131208539289e+0006L,
 435         1.720391071006963176836108026556547062980e+0007L,
 436         1.782614812922865190479394509487941920612e+0007L,
 437         9.708016389605273153536452032839879950155e+0006L,
 438         2.476495084688170096480215640962175140027e+0006L,
 439         2.363200660365585759668077790194604917187e+0005L,
 440         4.803239569848196077121203575704356936731e+0003L,
 441 };
 442 
 443 static const GENERIC pr6[13] = {        /* [1.28, 1.777...] */
 444         9.999999969777095495998606925524322559556e-0001L,
 445         5.825486719466194430503283824096872219216e+0001L,
 446         1.248155491637757281915184824965379905380e+0003L,
 447         1.302093199842358609321338417071710477615e+0004L,
 448         7.353835804186292782840961999810543016039e+0004L,
 449         2.356471661113686180549195092555751341757e+0005L,
 450         4.350553267429009581632987060942780847101e+0005L,
 451         4.588762661876600638719159826652389418235e+0005L,
 452         2.675796398548523436544221045225290128611e+0005L,
 453         8.077649557108971388298292919988449940464e+0004L,
 454         1.117640459221306873519068741664054573776e+0004L,
 455         5.544400072396814695175787511557757885585e+0002L,
 456         5.072550541191480498431289089905822910718e+0000L,
 457 };
 458 
 459 static const GENERIC ps6[13] = {        /* [1.28, 1.777...] */
 460         1.0e0L,
 461         5.832517925357165050639075848183613063291e+0001L,
 462         1.252144364743592128171256104364976466898e+0003L,
 463         1.310300234342216813579118022415585740772e+0004L,
 464         7.434667697093812197817292154032863632923e+0004L,
 465         2.398706595587719165726469002404004614711e+0005L,
 466         4.472737517625103157004869372427480602511e+0005L,
 467         4.786313523337761975294171429067037723611e+0005L,
 468         2.851161872872731228472536061865365370192e+0005L,
 469         8.891648269899148412331918021801385815586e+0004L,
 470         1.297097489535351517572978123584751042287e+0004L,
 471         7.096761640545975756202184143400469812618e+0002L,
 472         8.378049338590233325977702401733340820351e+0000L,
 473 };
 474 
 475 static const GENERIC sixteen = 16.0L;
 476 static const GENERIC huge = 1.0e30L;
 477 
 478 static GENERIC
 479 pzero(GENERIC x)
 480 {
 481         GENERIC s, r, t, z;
 482         int i;
 483 
 484         if (x > huge)
 485                 return (one);
 486 
 487         t = one / x;
 488         z = t * t;
 489 
 490         if (x > sixteen) {
 491                 r = z * pr0[11] + pr0[10];
 492                 s = ps0[10];
 493 
 494                 for (i = 9; i >= 0; i--) {
 495                         r = z * r + pr0[i];
 496                         s = z * s + ps0[i];
 497                 }
 498         } else if (x > eight) {
 499                 r = pr1[11];
 500                 s = ps1[11] + z * (ps1[12] + z * ps1[13]);
 501 
 502                 for (i = 10; i >= 0; i--) {
 503                         r = z * r + pr1[i];
 504                         s = z * s + ps1[i];
 505                 }
 506         } else if (x > five) {               /* x > 5.0 */
 507                 r = pr2[11];
 508                 s = ps2[11] + z * (ps2[12] + z * ps2[13]);
 509 
 510                 for (i = 10; i >= 0; i--) {
 511                         r = z * r + pr2[i];
 512                         s = z * s + ps2[i];
 513                 }
 514         } else if (x > 3.5L) {
 515                 r = pr3[12];
 516                 s = ps3[12];
 517 
 518                 for (i = 11; i >= 0; i--) {
 519                         r = z * r + pr3[i];
 520                         s = z * s + ps3[i];
 521                 }
 522         } else if (x > 2.5L) {
 523                 r = pr4[12];
 524                 s = ps4[12];
 525 
 526                 for (i = 11; i >= 0; i--) {
 527                         r = z * r + pr4[i];
 528                         s = z * s + ps4[i];
 529                 }
 530         } else if (x > (1.0L / 0.5625L)) {
 531                 r = pr5[12];
 532                 s = ps5[12];
 533 
 534                 for (i = 11; i >= 0; i--) {
 535                         r = z * r + pr5[i];
 536                         s = z * s + ps5[i];
 537                 }
 538         } else {                        /* assume x > 1.28 */
 539                 r = pr6[12];
 540                 s = ps6[12];
 541 
 542                 for (i = 11; i >= 0; i--) {
 543                         r = z * r + pr6[i];
 544                         s = z * s + ps6[i];
 545                 }
 546         }



 547 
 548         return (r / s);
 549 }
 550 static const GENERIC qr0[12] = {        /* [16, inf] */
 551         -1.249999999999999999999999999999999972972e-0001L,
 552         -1.425179595545670577414395762503991596897e+0002L,
 553         -6.312499645625970845534460257936222407219e+0004L,
 554         -1.411374326457208384315121243698814446848e+0007L,
 555         -1.735034212758873581410984757860787252842e+0009L,
 556         -1.199777647512789489421826342485055280680e+0011L,
 557         -4.596025334081655714499860409699100373644e+0012L,
 558         -9.262525628201284107792924477031653399187e+0013L,
 559         -8.858394728685039245344398842180662867639e+0014L,
 560         -3.267527953687534887623740622709505972113e+0015L,
 561         -2.664222971186311967587129347029450062019e+0015L,
 562         3.442464060723987869585180095344504100204e+0014L,
 563 };
 564 
 565 static const GENERIC qs0[11] = {
 566         1.0e0L,
 567         1.140729613936536461931516610003185687881e+0003L,
 568         5.056665510442299351009198186490085803580e+0005L,
 569         1.132041763825642787943941650522718199115e+0008L,
 570         1.394570111872581606392620678214246479767e+0010L,
 571         9.677945218152264789534431079563744378421e+0011L,
 572         3.731140327851536828225143058896348502096e+0013L,
 573         7.612785951064869291722846681020881676410e+0014L,
 574         7.476077016406764891730191004811863975940e+0015L,
 575         2.951246482613592035421503427100393831709e+0016L,
 576         3.108361803691811711136854587074302034901e+0016L,
 577 };
 578 
 579 static const GENERIC qr1[12] = {        /* [8, 16 ] */
 580         -1.249999999999999999997949010383433818157e-0001L,
 581         -9.051215166393822640636752244895124126934e+0001L,
 582         -2.620782703428148837671179031904208303947e+0004L,
 583         -3.975571261553504457766177974508785790884e+0006L,
 584         -3.479029330759311306270072218074074994090e+0008L,
 585         -1.823955008124268573036216746186239829089e+0010L,
 586         -5.765932697111801375765156029221568664435e+0011L,
 587         -1.079843680798742592954002192417934779114e+0013L,
 588         -1.146893630504592739082205764611581332897e+0014L,
 589         -6.367016059683898464936104447282880704182e+0014L,
 590         -1.583109041961213490464459111903484209098e+0015L,
 591         -1.230149555764242473103128650135795639412e+0015L,
 592 };
 593 
 594 static const GENERIC qs1[14] = {
 595         1.0e0L,
 596         7.246831508115058112438579847778014458432e+0002L,
 597         2.100854184439168518399383786306927037611e+0005L,
 598         3.192636418837951507430188285940994235122e+0007L,
 599         2.801558443383354674538443461124434216152e+0009L,
 600         1.475026997664373739293483927250653467487e+0011L,
 601         4.694486824913954608552363821799927145318e+0012L,
 602         8.890350100919200250838438709601547334021e+0013L,
 603         9.626844429082905144874701068760469752067e+0014L,
 604         5.541110744600460773528263862687521642140e+0015L,
 605         1.486500494789452556727470329232123096563e+0016L,
 606         1.415840104845959400365430773732093899210e+0016L,
 607         1.780866095241517418081312567239682336483e+0015L,
 608         -2.359230917384889357887631544079990129494e+0014L,
 609 };
 610 
 611 static const GENERIC qr2[12] = {        /* [5, 8] */
 612         -1.249999999999999531937744362527772181614e-0001L,
 613         -4.944373897356969774839375977239241573966e+0001L,
 614         -7.728449175433465285314261650078450473909e+0003L,
 615         -6.262574329612752346336901434651220705903e+0005L,
 616         -2.900948220220943306027235217424380672732e+0007L,
 617         -7.988719647634192770463917157562874119535e+0008L,
 618         -1.318228171927181389547760026626357012375e+0010L,
 619         -1.282439773983029245309263271945424928196e+0011L,
 620         -7.050925570827818040186149940257918845138e+0011L,
 621         -2.021751882573871990004205616874202684429e+0012L,
 622         -2.592939962400668552384333900573812635658e+0012L,
 623         -1.038267109518891262840601514932972850326e+0012L,
 624 };
 625 
 626 static const GENERIC qs2[14] = {
 627         1.0e0L,
 628         3.961358492885570003202784022894248952116e+0002L,
 629         6.205788738864701882828752634586510926968e+0004L,
 630         5.045715603932670286550673813011764406749e+0006L,
 631         2.349248611362658323353343389430968751429e+0008L,
 632         6.520244524415828635917683553721880063911e+0009L,
 633         1.089111211223507719337067159886281887722e+0011L,
 634         1.080406000905359867958779409414903018610e+0012L,
 635         6.135645280895514703514154680623769562148e+0012L,
 636         1.862433040246625874245867151368643668215e+0013L,
 637         2.667780805786648888840777888702193708994e+0013L,
 638         1.394401107289087774765300711809313112824e+0013L,
 639         1.093247500616320375562898297156722445484e+0012L,
 640         -7.228875530378928722826604216491493780775e+0010L,
 641 };
 642 
 643 static const GENERIC qr3[13] = {        /* [3.5 5] */
 644         -1.249999999999473067748420379578481661075e-0001L,
 645         -3.044549048635289351913574324803250977998e+0001L,
 646         -2.890081140649769078496693003524681440869e+0003L,
 647         -1.404922456817202235879343275330529107684e+0005L,
 648         -3.862746614385573443518177403617349281869e+0006L,
 649         -6.257517309110249049201133708911155047689e+0007L,
 650         -6.031451330920839916987079782727323477520e+0008L,
 651         -3.411542405173830611454025765755854382346e+0009L,
 652         -1.089392478149726672133014498723021526099e+0010L,
 653         -1.824934078420210941290140903415956782726e+0010L,
 654         -1.400780278304358710423481070486939531139e+0010L,
 655         -3.716484136064917363926635716743771092093e+0009L,
 656         -1.397591075296425529970434890954904331580e+0008L,
 657 };
 658 
 659 static const GENERIC qs3[13] = {
 660         1.0e0L,
 661         2.441498613904962049391000187014945858042e+0002L,
 662         2.326188882072370711500164222341514337043e+0004L,
 663         1.137138213121231338494977104659239578165e+0006L,
 664         3.152918070735662728722998452605364253517e+0007L,
 665         5.172877993426507259314270488444013595108e+0008L,
 666         5.083086439731669807455961078856470774115e+0009L,
 667         2.961842732066434123119325521139476909941e+0010L,
 668         9.912185866862440735829781856081353151390e+0010L,
 669         1.793560561251622234430564181567297983598e+0011L,
 670         1.577090119341228122525265108497940403073e+0011L,
 671         5.509910306780166194333889999985463681636e+0010L,
 672         4.761691134078874491202320181517936758141e+0009L,
 673 };
 674 
 675 static const GENERIC qr4[13] = {        /* [2.5 3.5] */
 676         -1.249999999928567734339745043490705340835e-0001L,
 677         -1.967201748731419063051601624435565528481e+0001L,
 678         -1.186329146714562236407099740615528170707e+0003L,
 679         -3.607736959222941810356301491152457934060e+0004L,
 680         -6.119200717978104904932828468575194267125e+0005L,
 681         -6.037847781158358226670305078652205586384e+0006L,
 682         -3.503558153336140359700536720393565984740e+0007L,
 683         -1.180196478268225718757218523746787309773e+0008L,
 684         -2.221860232085134915841426363505169680528e+0008L,
 685         -2.173372505452747585296176761701746236760e+0008L,
 686         -9.649364865061237558517730539506568013963e+0007L,
 687         -1.465429227847933034546039640094862650385e+0007L,
 688         -3.083003197920262085170581866246663380607e+0005L,
 689 };
 690 
 691 static const GENERIC qs4[13] = {        /* [2.5 3.5] */
 692         1.0e0L,
 693         1.579620773732259142752614142139986854055e+0002L,
 694         9.581372220329138733203879503753685054968e+0003L,
 695         2.939598672379108095776114131010825885308e+0005L,
 696         5.052183049314542218630341818692588448168e+0006L,
 697         5.083497695595206639433839326338971980149e+0007L,
 698         3.036385361800553388049719014005099206516e+0008L,
 699         1.067826481452753409910563785161661492137e+0009L,
 700         2.145644125557118044720741775125319669272e+0009L,
 701         2.324115615959719949363946673491552216799e+0009L,
 702         1.223262962112070757966959855619847011146e+0009L,
 703         2.569765553318495423738478585947110270709e+0008L,
 704         1.354744744299227127897905787732636565504e+0007L,
 705 };
 706 
 707 static const GENERIC qr5[13] = {        /* [1.777.., 2.5] */
 708         -1.249999995936639697637680428174576069971e-0001L,
 709         -1.260846055371311453485891923426489068315e+0001L,
 710         -4.772398467544467480801174330290141578895e+0002L,
 711         -8.939852599990298486613760833996490599724e+0003L,
 712         -9.184070787149542050979542226446134243197e+0004L,
 713         -5.406038945018274458362637897739280435171e+0005L,
 714         -1.845896544705190261018653728678171084418e+0006L,
 715         -3.613616990680809501878667570653308071547e+0006L,
 716         -3.908782978135693252252557720414348623779e+0006L,
 717         -2.173711022517323927109138170588442768176e+0006L,
 718         -5.431253130679918485836408549007856244495e+0005L,
 719         -4.591098546452684510082591587275940765959e+0004L,
 720         -5.244711364168207806835520057792229646578e+0002L,
 721 };
 722 
 723 static const GENERIC qs5[13] = {        /* [1.777.., 2.5] */
 724         1.0e0L,
 725         1.014536210851290878350892750972474861447e+0002L,
 726         3.875547510687135314064434160096139681076e+0003L,
 727         7.361913122670079814955259281995617732580e+0004L,
 728         7.720288944218771126581086539585529314636e+0005L,
 729         4.681529554446752496404431433608306558038e+0006L,
 730         1.667882621940503925455031252308367745820e+0007L,
 731         3.469403153761399881888272620855305156241e+0007L,
 732         4.096992047964210711867089384719947863019e+0007L,
 733         2.596804755829217449311530735959560630554e+0007L,
 734         7.983933774697889238154465064019410763845e+0006L,
 735         9.818133816979900819087242425280757938152e+0005L,
 736         3.061083930868694396013541535670745443560e+0004L,
 737 };
 738 
 739 static const GENERIC qr6[13] = {        /* [1.28, 1.777..] */
 740         -1.249999881577289001807137282824929082771e-0001L,
 741         -7.998273510053110759610810594119533619282e+0000L,
 742         -1.872481955335172543369089617771565632719e+0002L,
 743         -2.122116786726300805079874003303799646812e+0003L,
 744         -1.293850285839529282503178263484773478457e+0004L,
 745         -4.445024742266316181033354192262529356093e+0004L,
 746         -8.730161378334357767668344467356505347070e+0004L,
 747         -9.706222895172078442801444972505315054736e+0004L,
 748         -5.896325518259858270165531513618195321041e+0004L,
 749         -1.823172034368108822276420827074668832233e+0004L,
 750         -2.509304178635055926638833040337472387175e+0003L,
 751         -1.156608965715779237316769828941729964099e+0002L,
 752         -7.028005789650731396887346826397785210442e-0001L,
 753 };
 754 
 755 static const GENERIC qs6[13] = {        /* [1.28, 1.777..] */
 756         1.0e0L,
 757         6.457211085058064845601261321277721075900e+0001L,
 758         1.534005216588011210342824555136008682950e+0003L,
 759         1.777217999176441782593357660462379097171e+0004L,
 760         1.118372652642469468091084810263231199696e+0005L,
 761         4.015242433858461813142365748386473605294e+0005L,
 762         8.377081045517098645448616514388280497673e+0005L,
 763         1.011495020008010352575398009604164287337e+0006L,
 764         6.886722075290430568652227875200208955970e+0005L,
 765         2.504735189948021472047157148613171956537e+0005L,
 766         4.408138920171044846941001844352009817062e+0004L,
 767         3.105572178072115145673058722853640854884e+0003L,
 768         5.588294821118916113437396504573817033678e+0001L,
 769 };
 770 
 771 static GENERIC
 772 qzero(GENERIC x)
 773 {
 774         GENERIC s, r, t, z;
 775         int i;
 776 
 777         if (x > huge)
 778                 return (-0.125L / x);
 779 
 780         t = one / x;
 781         z = t * t;
 782 
 783         if (x > sixteen) {
 784                 r = z * qr0[11] + qr0[10];
 785                 s = qs0[10];
 786 
 787                 for (i = 9; i >= 0; i--) {
 788                         r = z * r + qr0[i];
 789                         s = z * s + qs0[i];
 790                 }
 791         } else if (x > eight) {
 792                 r = qr1[11];
 793                 s = qs1[11] + z * (qs1[12] + z * qs1[13]);
 794 
 795                 for (i = 10; i >= 0; i--) {
 796                         r = z * r + qr1[i];
 797                         s = z * s + qs1[i];
 798                 }
 799         } else if (x > five) {       /* assume x > 5.0 */
 800                 r = qr2[11];
 801                 s = qs2[11] + z * (qs2[12] + z * qs2[13]);
 802 
 803                 for (i = 10; i >= 0; i--) {
 804                         r = z * r + qr2[i];
 805                         s = z * s + qs2[i];
 806                 }
 807         } else if (x > 3.5L) {
 808                 r = qr3[12];
 809                 s = qs3[12];
 810 
 811                 for (i = 11; i >= 0; i--) {
 812                         r = z * r + qr3[i];
 813                         s = z * s + qs3[i];
 814                 }
 815         } else if (x > 2.5L) {
 816                 r = qr4[12];
 817                 s = qs4[12];
 818 
 819                 for (i = 11; i >= 0; i--) {
 820                         r = z * r + qr4[i];
 821                         s = z * s + qs4[i];
 822                 }
 823         } else if (x > (1.0L / 0.5625L)) {
 824                 r = qr5[12];
 825                 s = qs5[12];
 826 
 827                 for (i = 11; i >= 0; i--) {
 828                         r = z * r + qr5[i];
 829                         s = z * s + qs5[i];
 830                 }
 831         } else {                        /* assume x > 1.28 */
 832                 r = qr6[12];
 833                 s = qs6[12];
 834 
 835                 for (i = 11; i >= 0; i--) {
 836                         r = z * r + qr6[i];
 837                         s = z * s + qs6[i];
 838                 }
 839         }
 840 
 841         return (t * (r / s));
 842 }