sysdeps/ieee754/ldbl-96/e_j1l.c

   1 /*
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10  */
  11
  12 /* Long double expansions are
  13   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
  14   and are incorporated herein by permission of the author.  The author
  15   reserves the right to distribute this material elsewhere under different
  16   copying permissions.  These modifications are distributed here under
  17   the following terms:
  18
  19     This library is free software; you can redistribute it and/or
  20     modify it under the terms of the GNU Lesser General Public
  21     License as published by the Free Software Foundation; either
  22     version 2.1 of the License, or (at your option) any later version.
  23
  24     This library is distributed in the hope that it will be useful,
  25     but WITHOUT ANY WARRANTY; without even the implied warranty of
  26     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  27     Lesser General Public License for more details.
  28
  29     You should have received a copy of the GNU Lesser General Public
  30     License along with this library; if not, see
  31     <http://www.gnu.org/licenses/>.  */
  32
  33 /* __ieee754_j1(x), __ieee754_y1(x)
  34  * Bessel function of the first and second kinds of order zero.
  35  * Method -- j1(x):
  36  *      1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
  37  *      2. Reduce x to |x| since j1(x)=-j1(-x),  and
  38  *         for x in (0,2)
  39  *              j1(x) = x/2 + x*z*R0/S0,  where z = x*x;
  40  *         for x in (2,inf)
  41  *              j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
  42  *              y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
  43  *         where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
  44  *         as follow:
  45  *              cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
  46  *                      =  1/sqrt(2) * (sin(x) - cos(x))
  47  *              sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
  48  *                      = -1/sqrt(2) * (sin(x) + cos(x))
  49  *         (To avoid cancellation, use
  50  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
  51  *          to compute the worse one.)
  52  *
  53  *      3 Special cases
  54  *              j1(nan)= nan
  55  *              j1(0) = 0
  56  *              j1(inf) = 0
  57  *
  58  * Method -- y1(x):
  59  *      1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
  60  *      2. For x<2.
  61  *         Since
  62  *              y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
  63  *         therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
  64  *         We use the following function to approximate y1,
  65  *              y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
  66  *         Note: For tiny x, 1/x dominate y1 and hence
  67  *              y1(tiny) = -2/pi/tiny
  68  *      3. For x>=2.
  69  *              y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
  70  *         where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
  71  *         by method mentioned above.
  72  */
  73
  74 #include <errno.h>
  75 #include <float.h>
  76 #include <math.h>
  77 #include <math_private.h>
  78
  79 static long double pone (long double), qone (long double);
  80
  81 static const long double
  82   huge = 1e4930L,
  83  one = 1.0L,
  84  invsqrtpi = 5.6418958354775628694807945156077258584405e-1L,
  85   tpi =  6.3661977236758134307553505349005744813784e-1L,
  86
  87   /* J1(x) = .5 x + x x^2 R(x^2) / S(x^2)
  88      0 <= x <= 2
  89      Peak relative error 4.5e-21 */
  90 R[5] = {
  91     -9.647406112428107954753770469290757756814E7L,
  92     2.686288565865230690166454005558203955564E6L,
  93     -3.689682683905671185891885948692283776081E4L,
  94     2.195031194229176602851429567792676658146E2L,
  95     -5.124499848728030297902028238597308971319E-1L,
  96 },
  97
  98   S[4] =
  99 {
 100   1.543584977988497274437410333029029035089E9L,
 101   2.133542369567701244002565983150952549520E7L,
 102   1.394077011298227346483732156167414670520E5L,
 103   5.252401789085732428842871556112108446506E2L,
 104   /*  1.000000000000000000000000000000000000000E0L, */
 105 };
 106
 107 static const long double zero = 0.0;
 108
 109
 110 long double
 111 __ieee754_j1l (long double x)
 112 {
 113   long double z, c, r, s, ss, cc, u, v, y;
 114   int32_t ix;
 115   uint32_t se;
 116
 117   GET_LDOUBLE_EXP (se, x);
 118   ix = se & 0x7fff;
 119   if (__glibc_unlikely (ix >= 0x7fff))
 120     return one / x;
 121   y = fabsl (x);
 122   if (ix >= 0x4000)
 123     {                           /* |x| >= 2.0 */
 124       __sincosl (y, &s, &c);
 125       ss = -s - c;
 126       cc = s - c;
 127       if (ix < 0x7ffe)
 128         {                       /* make sure y+y not overflow */
 129           z = __cosl (y + y);
 130           if ((s * c) > zero)
 131             cc = z / ss;
 132           else
 133             ss = z / cc;
 134         }
 135       /*
 136        * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
 137        * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
 138        */
 139       if (__glibc_unlikely (ix > 0x4080))
 140         z = (invsqrtpi * cc) / sqrtl (y);
 141       else
 142         {
 143           u = pone (y);
 144           v = qone (y);
 145           z = invsqrtpi * (u * cc - v * ss) / sqrtl (y);
 146         }
 147       if (se & 0x8000)
 148         return -z;
 149       else
 150         return z;
 151     }
 152   if (__glibc_unlikely (ix < 0x3fde))       /* |x| < 2^-33 */
 153     {
 154       if (huge + x > one)               /* inexact if x!=0 necessary */
 155         {
 156           long double ret = 0.5 * x;
 157           math_check_force_underflow (ret);
 158           if (ret == 0 && x != 0)
 159             __set_errno (ERANGE);
 160           return ret;
 161         }
 162     }
 163   z = x * x;
 164   r = z * (R[0] + z * (R[1]+ z * (R[2] + z * (R[3] + z * R[4]))));
 165   s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z)));
 166   r *= x;
 167   return (x * 0.5 + r / s);
 168 }
 169 strong_alias (__ieee754_j1l, __j1l_finite)
 170
 171
 172 /* Y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + x R(x^2)
 173    0 <= x <= 2
 174    Peak relative error 2.3e-23 */
 175 static const long double U0[6] = {
 176   -5.908077186259914699178903164682444848615E10L,
 177   1.546219327181478013495975514375773435962E10L,
 178   -6.438303331169223128870035584107053228235E8L,
 179   9.708540045657182600665968063824819371216E6L,
 180   -6.138043997084355564619377183564196265471E4L,
 181   1.418503228220927321096904291501161800215E2L,
 182 };
 183 static const long double V0[5] = {
 184   3.013447341682896694781964795373783679861E11L,
 185   4.669546565705981649470005402243136124523E9L,
 186   3.595056091631351184676890179233695857260E7L,
 187   1.761554028569108722903944659933744317994E5L,
 188   5.668480419646516568875555062047234534863E2L,
 189   /*  1.000000000000000000000000000000000000000E0L, */
 190 };
 191
 192
 193 long double
 194 __ieee754_y1l (long double x)
 195 {
 196   long double z, s, c, ss, cc, u, v;
 197   int32_t ix;
 198   uint32_t se, i0, i1;
 199
 200   GET_LDOUBLE_WORDS (se, i0, i1, x);
 201   ix = se & 0x7fff;
 202   /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
 203   if (__glibc_unlikely (se & 0x8000))
 204     return zero / (zero * x);
 205   if (__glibc_unlikely (ix >= 0x7fff))
 206     return one / (x + x * x);
 207   if (__glibc_unlikely ((i0 | i1) == 0))
 208     return -HUGE_VALL + x;  /* -inf and overflow exception.  */
 209   if (ix >= 0x4000)
 210     {                           /* |x| >= 2.0 */
 211       __sincosl (x, &s, &c);
 212       ss = -s - c;
 213       cc = s - c;
 214       if (ix < 0x7ffe)
 215         {                       /* make sure x+x not overflow */
 216           z = __cosl (x + x);
 217           if ((s * c) > zero)
 218             cc = z / ss;
 219           else
 220             ss = z / cc;
 221         }
 222       /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
 223        * where x0 = x-3pi/4
 224        *      Better formula:
 225        *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
 226        *                      =  1/sqrt(2) * (sin(x) - cos(x))
 227        *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 228        *                      = -1/sqrt(2) * (cos(x) + sin(x))
 229        * To avoid cancellation, use
 230        *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 231        * to compute the worse one.
 232        */
 233       if (__glibc_unlikely (ix > 0x4080))
 234         z = (invsqrtpi * ss) / sqrtl (x);
 235       else
 236         {
 237           u = pone (x);
 238           v = qone (x);
 239           z = invsqrtpi * (u * ss + v * cc) / sqrtl (x);
 240         }
 241       return z;
 242     }
 243   if (__glibc_unlikely (ix <= 0x3fbe))
 244     {                           /* x < 2**-65 */
 245       z = -tpi / x;
 246       if (isinf (z))
 247         __set_errno (ERANGE);
 248       return z;
 249     }
 250   z = x * x;
 251  u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * (U0[4] + z * U0[5]))));
 252  v = V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * (V0[4] + z))));
 253   return (x * (u / v) +
 254           tpi * (__ieee754_j1l (x) * __ieee754_logl (x) - one / x));
 255 }
 256 strong_alias (__ieee754_y1l, __y1l_finite)
 257
 258
 259 /* For x >= 8, the asymptotic expansions of pone is
 260  *      1 + 15/128 s^2 - 4725/2^15 s^4 - ...,   where s = 1/x.
 261  * We approximate pone by
 262  *      pone(x) = 1 + (R/S)
 263  */
 264
 265 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 266    P1(x) = 1 + z^2 R(z^2), z=1/x
 267    8 <= x <= inf  (0 <= z <= 0.125)
 268    Peak relative error 5.2e-22  */
 269
 270 static const long double pr8[7] = {
 271   8.402048819032978959298664869941375143163E-9L,
 272   1.813743245316438056192649247507255996036E-6L,
 273   1.260704554112906152344932388588243836276E-4L,
 274   3.439294839869103014614229832700986965110E-3L,
 275   3.576910849712074184504430254290179501209E-2L,
 276   1.131111483254318243139953003461511308672E-1L,
 277   4.480715825681029711521286449131671880953E-2L,
 278 };
 279 static const long double ps8[6] = {
 280   7.169748325574809484893888315707824924354E-8L,
 281   1.556549720596672576431813934184403614817E-5L,
 282   1.094540125521337139209062035774174565882E-3L,
 283   3.060978962596642798560894375281428805840E-2L,
 284   3.374146536087205506032643098619414507024E-1L,
 285   1.253830208588979001991901126393231302559E0L,
 286   /* 1.000000000000000000000000000000000000000E0L, */
 287 };
 288
 289 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 290    P1(x) = 1 + z^2 R(z^2), z=1/x
 291    4.54541015625 <= x <= 8
 292    Peak relative error 7.7e-22  */
 293 static const long double pr5[7] = {
 294   4.318486887948814529950980396300969247900E-7L,
 295   4.715341880798817230333360497524173929315E-5L,
 296   1.642719430496086618401091544113220340094E-3L,
 297   2.228688005300803935928733750456396149104E-2L,
 298   1.142773760804150921573259605730018327162E-1L,
 299   1.755576530055079253910829652698703791957E-1L,
 300   3.218803858282095929559165965353784980613E-2L,
 301 };
 302 static const long double ps5[6] = {
 303   3.685108812227721334719884358034713967557E-6L,
 304   4.069102509511177498808856515005792027639E-4L,
 305   1.449728676496155025507893322405597039816E-2L,
 306   2.058869213229520086582695850441194363103E-1L,
 307   1.164890985918737148968424972072751066553E0L,
 308   2.274776933457009446573027260373361586841E0L,
 309   /*  1.000000000000000000000000000000000000000E0L,*/
 310 };
 311
 312 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 313    P1(x) = 1 + z^2 R(z^2), z=1/x
 314    2.85711669921875 <= x <= 4.54541015625
 315    Peak relative error 6.5e-21  */
 316 static const long double pr3[7] = {
 317   1.265251153957366716825382654273326407972E-5L,
 318   8.031057269201324914127680782288352574567E-4L,
 319   1.581648121115028333661412169396282881035E-2L,
 320   1.179534658087796321928362981518645033967E-1L,
 321   3.227936912780465219246440724502790727866E-1L,
 322   2.559223765418386621748404398017602935764E-1L,
 323   2.277136933287817911091370397134882441046E-2L,
 324 };
 325 static const long double ps3[6] = {
 326   1.079681071833391818661952793568345057548E-4L,
 327   6.986017817100477138417481463810841529026E-3L,
 328   1.429403701146942509913198539100230540503E-1L,
 329   1.148392024337075609460312658938700765074E0L,
 330   3.643663015091248720208251490291968840882E0L,
 331   3.990702269032018282145100741746633960737E0L,
 332   /* 1.000000000000000000000000000000000000000E0L, */
 333 };
 334
 335 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 336    P1(x) = 1 + z^2 R(z^2), z=1/x
 337    2 <= x <= 2.85711669921875
 338    Peak relative error 3.5e-21  */
 339 static const long double pr2[7] = {
 340   2.795623248568412225239401141338714516445E-4L,
 341   1.092578168441856711925254839815430061135E-2L,
 342   1.278024620468953761154963591853679640560E-1L,
 343   5.469680473691500673112904286228351988583E-1L,
 344   8.313769490922351300461498619045639016059E-1L,
 345   3.544176317308370086415403567097130611468E-1L,
 346   1.604142674802373041247957048801599740644E-2L,
 347 };
 348 static const long double ps2[6] = {
 349   2.385605161555183386205027000675875235980E-3L,
 350   9.616778294482695283928617708206967248579E-2L,
 351   1.195215570959693572089824415393951258510E0L,
 352   5.718412857897054829999458736064922974662E0L,
 353   1.065626298505499086386584642761602177568E1L,
 354   6.809140730053382188468983548092322151791E0L,
 355  /* 1.000000000000000000000000000000000000000E0L, */
 356 };
 357
 358
 359 static long double
 360 pone (long double x)
 361 {
 362   const long double *p, *q;
 363   long double z, r, s;
 364   int32_t ix;
 365   uint32_t se, i0, i1;
 366
 367   GET_LDOUBLE_WORDS (se, i0, i1, x);
 368   ix = se & 0x7fff;
 369   /* ix >= 0x4000 for all calls to this function.  */
 370   if (ix >= 0x4002) /* x >= 8 */
 371     {
 372       p = pr8;
 373       q = ps8;
 374     }
 375   else
 376     {
 377       i1 = (ix << 16) | (i0 >> 16);
 378       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 379         {
 380           p = pr5;
 381           q = ps5;
 382         }
 383       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 384         {
 385           p = pr3;
 386           q = ps3;
 387         }
 388       else      /* x >= 2 */
 389         {
 390           p = pr2;
 391           q = ps2;
 392         }
 393     }
 394   z = one / (x * x);
 395  r = p[0] + z * (p[1] +
 396                  z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 397  s = q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
 398   return one + z * r / s;
 399 }
 400
 401
 402 /* For x >= 8, the asymptotic expansions of qone is
 403  *      3/8 s - 105/1024 s^3 - ..., where s = 1/x.
 404  * We approximate pone by
 405  *      qone(x) = s*(0.375 + (R/S))
 406  */
 407
 408 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 409    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 410    8 <= x <= inf
 411    Peak relative error 8.3e-22 */
 412
 413 static const long double qr8[7] = {
 414   -5.691925079044209246015366919809404457380E-10L,
 415   -1.632587664706999307871963065396218379137E-7L,
 416   -1.577424682764651970003637263552027114600E-5L,
 417   -6.377627959241053914770158336842725291713E-4L,
 418   -1.087408516779972735197277149494929568768E-2L,
 419   -6.854943629378084419631926076882330494217E-2L,
 420   -1.055448290469180032312893377152490183203E-1L,
 421 };
 422 static const long double qs8[7] = {
 423   5.550982172325019811119223916998393907513E-9L,
 424   1.607188366646736068460131091130644192244E-6L,
 425   1.580792530091386496626494138334505893599E-4L,
 426   6.617859900815747303032860443855006056595E-3L,
 427   1.212840547336984859952597488863037659161E-1L,
 428   9.017885953937234900458186716154005541075E-1L,
 429   2.201114489712243262000939120146436167178E0L,
 430   /* 1.000000000000000000000000000000000000000E0L, */
 431 };
 432
 433 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 434    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 435    4.54541015625 <= x <= 8
 436    Peak relative error 4.1e-22 */
 437 static const long double qr5[7] = {
 438   -6.719134139179190546324213696633564965983E-8L,
 439   -9.467871458774950479909851595678622044140E-6L,
 440   -4.429341875348286176950914275723051452838E-4L,
 441   -8.539898021757342531563866270278505014487E-3L,
 442   -6.818691805848737010422337101409276287170E-2L,
 443   -1.964432669771684034858848142418228214855E-1L,
 444   -1.333896496989238600119596538299938520726E-1L,
 445 };
 446 static const long double qs5[7] = {
 447   6.552755584474634766937589285426911075101E-7L,
 448   9.410814032118155978663509073200494000589E-5L,
 449   4.561677087286518359461609153655021253238E-3L,
 450   9.397742096177905170800336715661091535805E-2L,
 451   8.518538116671013902180962914473967738771E-1L,
 452   3.177729183645800174212539541058292579009E0L,
 453   4.006745668510308096259753538973038902990E0L,
 454   /* 1.000000000000000000000000000000000000000E0L, */
 455 };
 456
 457 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 458    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 459    2.85711669921875 <= x <= 4.54541015625
 460    Peak relative error 2.2e-21 */
 461 static const long double qr3[7] = {
 462   -3.618746299358445926506719188614570588404E-6L,
 463   -2.951146018465419674063882650970344502798E-4L,
 464   -7.728518171262562194043409753656506795258E-3L,
 465   -8.058010968753999435006488158237984014883E-2L,
 466   -3.356232856677966691703904770937143483472E-1L,
 467   -4.858192581793118040782557808823460276452E-1L,
 468   -1.592399251246473643510898335746432479373E-1L,
 469 };
 470 static const long double qs3[7] = {
 471   3.529139957987837084554591421329876744262E-5L,
 472   2.973602667215766676998703687065066180115E-3L,
 473   8.273534546240864308494062287908662592100E-2L,
 474   9.613359842126507198241321110649974032726E-1L,
 475   4.853923697093974370118387947065402707519E0L,
 476   1.002671608961669247462020977417828796933E1L,
 477   7.028927383922483728931327850683151410267E0L,
 478   /* 1.000000000000000000000000000000000000000E0L, */
 479 };
 480
 481 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 482    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 483    2 <= x <= 2.85711669921875
 484    Peak relative error 6.9e-22 */
 485 static const long double qr2[7] = {
 486   -1.372751603025230017220666013816502528318E-4L,
 487   -6.879190253347766576229143006767218972834E-3L,
 488   -1.061253572090925414598304855316280077828E-1L,
 489   -6.262164224345471241219408329354943337214E-1L,
 490   -1.423149636514768476376254324731437473915E0L,
 491   -1.087955310491078933531734062917489870754E0L,
 492   -1.826821119773182847861406108689273719137E-1L,
 493 };
 494 static const long double qs2[7] = {
 495   1.338768933634451601814048220627185324007E-3L,
 496   7.071099998918497559736318523932241901810E-2L,
 497   1.200511429784048632105295629933382142221E0L,
 498   8.327301713640367079030141077172031825276E0L,
 499   2.468479301872299311658145549931764426840E1L,
 500   2.961179686096262083509383820557051621644E1L,
 501   1.201402313144305153005639494661767354977E1L,
 502  /* 1.000000000000000000000000000000000000000E0L, */
 503 };
 504
 505
 506 static long double
 507 qone (long double x)
 508 {
 509   const long double *p, *q;
 510   long double s, r, z;
 511   int32_t ix;
 512   uint32_t se, i0, i1;
 513
 514   GET_LDOUBLE_WORDS (se, i0, i1, x);
 515   ix = se & 0x7fff;
 516   /* ix >= 0x4000 for all calls to this function.  */
 517   if (ix >= 0x4002)             /* x >= 8 */
 518     {
 519       p = qr8;
 520       q = qs8;
 521     }
 522   else
 523     {
 524       i1 = (ix << 16) | (i0 >> 16);
 525       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 526         {
 527           p = qr5;
 528           q = qs5;
 529         }
 530       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 531         {
 532           p = qr3;
 533           q = qs3;
 534         }
 535       else      /* x >= 2 */
 536         {
 537           p = qr2;
 538           q = qs2;
 539         }
 540     }
 541   z = one / (x * x);
 542   r =
 543     p[0] + z * (p[1] +
 544                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 545   s =
 546     q[0] + z * (q[1] +
 547                 z * (q[2] +
 548                      z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
 549   return (.375 + z * r / s) / x;
 550 }