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     <https://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 #include <math-underflow.h>
  79 #include <libm-alias-finite.h>
  80
  81 static long double pone (long double), qone (long double);
  82
  83 static const long double
  84   huge = 1e4930L,
  85  one = 1.0L,
  86  invsqrtpi = 5.6418958354775628694807945156077258584405e-1L,
  87   tpi =  6.3661977236758134307553505349005744813784e-1L,
  88
  89   /* J1(x) = .5 x + x x^2 R(x^2) / S(x^2)
  90      0 <= x <= 2
  91      Peak relative error 4.5e-21 */
  92 R[5] = {
  93     -9.647406112428107954753770469290757756814E7L,
  94     2.686288565865230690166454005558203955564E6L,
  95     -3.689682683905671185891885948692283776081E4L,
  96     2.195031194229176602851429567792676658146E2L,
  97     -5.124499848728030297902028238597308971319E-1L,
  98 },
  99
 100   S[4] =
 101 {
 102   1.543584977988497274437410333029029035089E9L,
 103   2.133542369567701244002565983150952549520E7L,
 104   1.394077011298227346483732156167414670520E5L,
 105   5.252401789085732428842871556112108446506E2L,
 106   /*  1.000000000000000000000000000000000000000E0L, */
 107 };
 108
 109 static const long double zero = 0.0;
 110
 111
 112 long double
 113 __ieee754_j1l (long double x)
 114 {
 115   long double z, c, r, s, ss, cc, u, v, y;
 116   int32_t ix;
 117   uint32_t se;
 118
 119   GET_LDOUBLE_EXP (se, x);
 120   ix = se & 0x7fff;
 121   if (__glibc_unlikely (ix >= 0x7fff))
 122     return one / x;
 123   y = fabsl (x);
 124   if (ix >= 0x4000)
 125     {                           /* |x| >= 2.0 */
 126       __sincosl (y, &s, &c);
 127       ss = -s - c;
 128       cc = s - c;
 129       if (ix < 0x7ffe)
 130         {                       /* make sure y+y not overflow */
 131           z = __cosl (y + y);
 132           if ((s * c) > zero)
 133             cc = z / ss;
 134           else
 135             ss = z / cc;
 136         }
 137       /*
 138        * j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
 139        * y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
 140        */
 141       if (__glibc_unlikely (ix > 0x408e))
 142         z = (invsqrtpi * cc) / sqrtl (y);
 143       else
 144         {
 145           u = pone (y);
 146           v = qone (y);
 147           z = invsqrtpi * (u * cc - v * ss) / sqrtl (y);
 148         }
 149       if (se & 0x8000)
 150         return -z;
 151       else
 152         return z;
 153     }
 154   if (__glibc_unlikely (ix < 0x3fde))       /* |x| < 2^-33 */
 155     {
 156       if (huge + x > one)               /* inexact if x!=0 necessary */
 157         {
 158           long double ret = 0.5 * x;
 159           math_check_force_underflow (ret);
 160           if (ret == 0 && x != 0)
 161             __set_errno (ERANGE);
 162           return ret;
 163         }
 164     }
 165   z = x * x;
 166   r = z * (R[0] + z * (R[1]+ z * (R[2] + z * (R[3] + z * R[4]))));
 167   s = S[0] + z * (S[1] + z * (S[2] + z * (S[3] + z)));
 168   r *= x;
 169   return (x * 0.5 + r / s);
 170 }
 171 libm_alias_finite (__ieee754_j1l, __j1l)
 172
 173
 174 /* Y1(x) = 2/pi * (log(x) * j1(x) - 1/x) + x R(x^2)
 175    0 <= x <= 2
 176    Peak relative error 2.3e-23 */
 177 static const long double U0[6] = {
 178   -5.908077186259914699178903164682444848615E10L,
 179   1.546219327181478013495975514375773435962E10L,
 180   -6.438303331169223128870035584107053228235E8L,
 181   9.708540045657182600665968063824819371216E6L,
 182   -6.138043997084355564619377183564196265471E4L,
 183   1.418503228220927321096904291501161800215E2L,
 184 };
 185 static const long double V0[5] = {
 186   3.013447341682896694781964795373783679861E11L,
 187   4.669546565705981649470005402243136124523E9L,
 188   3.595056091631351184676890179233695857260E7L,
 189   1.761554028569108722903944659933744317994E5L,
 190   5.668480419646516568875555062047234534863E2L,
 191   /*  1.000000000000000000000000000000000000000E0L, */
 192 };
 193
 194
 195 long double
 196 __ieee754_y1l (long double x)
 197 {
 198   long double z, s, c, ss, cc, u, v;
 199   int32_t ix;
 200   uint32_t se, i0, i1;
 201
 202   GET_LDOUBLE_WORDS (se, i0, i1, x);
 203   ix = se & 0x7fff;
 204   /* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
 205   if (__glibc_unlikely (se & 0x8000))
 206     return zero / (zero * x);
 207   if (__glibc_unlikely (ix >= 0x7fff))
 208     return one / (x + x * x);
 209   if (__glibc_unlikely ((i0 | i1) == 0))
 210     return -HUGE_VALL + x;  /* -inf and overflow exception.  */
 211   if (ix >= 0x4000)
 212     {                           /* |x| >= 2.0 */
 213       __sincosl (x, &s, &c);
 214       ss = -s - c;
 215       cc = s - c;
 216       if (ix < 0x7ffe)
 217         {                       /* make sure x+x not overflow */
 218           z = __cosl (x + x);
 219           if ((s * c) > zero)
 220             cc = z / ss;
 221           else
 222             ss = z / cc;
 223         }
 224       /* y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x0)+q1(x)*cos(x0))
 225        * where x0 = x-3pi/4
 226        *      Better formula:
 227        *              cos(x0) = cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
 228        *                      =  1/sqrt(2) * (sin(x) - cos(x))
 229        *              sin(x0) = sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
 230        *                      = -1/sqrt(2) * (cos(x) + sin(x))
 231        * To avoid cancellation, use
 232        *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 233        * to compute the worse one.
 234        */
 235       if (__glibc_unlikely (ix > 0x408e))
 236         z = (invsqrtpi * ss) / sqrtl (x);
 237       else
 238         {
 239           u = pone (x);
 240           v = qone (x);
 241           z = invsqrtpi * (u * ss + v * cc) / sqrtl (x);
 242         }
 243       return z;
 244     }
 245   if (__glibc_unlikely (ix <= 0x3fbe))
 246     {                           /* x < 2**-65 */
 247       z = -tpi / x;
 248       if (isinf (z))
 249         __set_errno (ERANGE);
 250       return z;
 251     }
 252   z = x * x;
 253  u = U0[0] + z * (U0[1] + z * (U0[2] + z * (U0[3] + z * (U0[4] + z * U0[5]))));
 254  v = V0[0] + z * (V0[1] + z * (V0[2] + z * (V0[3] + z * (V0[4] + z))));
 255   return (x * (u / v) +
 256           tpi * (__ieee754_j1l (x) * __ieee754_logl (x) - one / x));
 257 }
 258 libm_alias_finite (__ieee754_y1l, __y1l)
 259
 260
 261 /* For x >= 8, the asymptotic expansions of pone is
 262  *      1 + 15/128 s^2 - 4725/2^15 s^4 - ...,   where s = 1/x.
 263  * We approximate pone by
 264  *      pone(x) = 1 + (R/S)
 265  */
 266
 267 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 268    P1(x) = 1 + z^2 R(z^2), z=1/x
 269    8 <= x <= inf  (0 <= z <= 0.125)
 270    Peak relative error 5.2e-22  */
 271
 272 static const long double pr8[7] = {
 273   8.402048819032978959298664869941375143163E-9L,
 274   1.813743245316438056192649247507255996036E-6L,
 275   1.260704554112906152344932388588243836276E-4L,
 276   3.439294839869103014614229832700986965110E-3L,
 277   3.576910849712074184504430254290179501209E-2L,
 278   1.131111483254318243139953003461511308672E-1L,
 279   4.480715825681029711521286449131671880953E-2L,
 280 };
 281 static const long double ps8[6] = {
 282   7.169748325574809484893888315707824924354E-8L,
 283   1.556549720596672576431813934184403614817E-5L,
 284   1.094540125521337139209062035774174565882E-3L,
 285   3.060978962596642798560894375281428805840E-2L,
 286   3.374146536087205506032643098619414507024E-1L,
 287   1.253830208588979001991901126393231302559E0L,
 288   /* 1.000000000000000000000000000000000000000E0L, */
 289 };
 290
 291 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 292    P1(x) = 1 + z^2 R(z^2), z=1/x
 293    4.54541015625 <= x <= 8
 294    Peak relative error 7.7e-22  */
 295 static const long double pr5[7] = {
 296   4.318486887948814529950980396300969247900E-7L,
 297   4.715341880798817230333360497524173929315E-5L,
 298   1.642719430496086618401091544113220340094E-3L,
 299   2.228688005300803935928733750456396149104E-2L,
 300   1.142773760804150921573259605730018327162E-1L,
 301   1.755576530055079253910829652698703791957E-1L,
 302   3.218803858282095929559165965353784980613E-2L,
 303 };
 304 static const long double ps5[6] = {
 305   3.685108812227721334719884358034713967557E-6L,
 306   4.069102509511177498808856515005792027639E-4L,
 307   1.449728676496155025507893322405597039816E-2L,
 308   2.058869213229520086582695850441194363103E-1L,
 309   1.164890985918737148968424972072751066553E0L,
 310   2.274776933457009446573027260373361586841E0L,
 311   /*  1.000000000000000000000000000000000000000E0L,*/
 312 };
 313
 314 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 315    P1(x) = 1 + z^2 R(z^2), z=1/x
 316    2.85711669921875 <= x <= 4.54541015625
 317    Peak relative error 6.5e-21  */
 318 static const long double pr3[7] = {
 319   1.265251153957366716825382654273326407972E-5L,
 320   8.031057269201324914127680782288352574567E-4L,
 321   1.581648121115028333661412169396282881035E-2L,
 322   1.179534658087796321928362981518645033967E-1L,
 323   3.227936912780465219246440724502790727866E-1L,
 324   2.559223765418386621748404398017602935764E-1L,
 325   2.277136933287817911091370397134882441046E-2L,
 326 };
 327 static const long double ps3[6] = {
 328   1.079681071833391818661952793568345057548E-4L,
 329   6.986017817100477138417481463810841529026E-3L,
 330   1.429403701146942509913198539100230540503E-1L,
 331   1.148392024337075609460312658938700765074E0L,
 332   3.643663015091248720208251490291968840882E0L,
 333   3.990702269032018282145100741746633960737E0L,
 334   /* 1.000000000000000000000000000000000000000E0L, */
 335 };
 336
 337 /* J1(x) cosX + Y1(x) sinX  =  sqrt( 2/(pi x)) P1(x)
 338    P1(x) = 1 + z^2 R(z^2), z=1/x
 339    2 <= x <= 2.85711669921875
 340    Peak relative error 3.5e-21  */
 341 static const long double pr2[7] = {
 342   2.795623248568412225239401141338714516445E-4L,
 343   1.092578168441856711925254839815430061135E-2L,
 344   1.278024620468953761154963591853679640560E-1L,
 345   5.469680473691500673112904286228351988583E-1L,
 346   8.313769490922351300461498619045639016059E-1L,
 347   3.544176317308370086415403567097130611468E-1L,
 348   1.604142674802373041247957048801599740644E-2L,
 349 };
 350 static const long double ps2[6] = {
 351   2.385605161555183386205027000675875235980E-3L,
 352   9.616778294482695283928617708206967248579E-2L,
 353   1.195215570959693572089824415393951258510E0L,
 354   5.718412857897054829999458736064922974662E0L,
 355   1.065626298505499086386584642761602177568E1L,
 356   6.809140730053382188468983548092322151791E0L,
 357  /* 1.000000000000000000000000000000000000000E0L, */
 358 };
 359
 360
 361 static long double
 362 pone (long double x)
 363 {
 364   const long double *p, *q;
 365   long double z, r, s;
 366   int32_t ix;
 367   uint32_t se, i0, i1;
 368
 369   GET_LDOUBLE_WORDS (se, i0, i1, x);
 370   ix = se & 0x7fff;
 371   /* ix >= 0x4000 for all calls to this function.  */
 372   if (ix >= 0x4002) /* x >= 8 */
 373     {
 374       p = pr8;
 375       q = ps8;
 376     }
 377   else
 378     {
 379       i1 = (ix << 16) | (i0 >> 16);
 380       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 381         {
 382           p = pr5;
 383           q = ps5;
 384         }
 385       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 386         {
 387           p = pr3;
 388           q = ps3;
 389         }
 390       else      /* x >= 2 */
 391         {
 392           p = pr2;
 393           q = ps2;
 394         }
 395     }
 396   z = one / (x * x);
 397  r = p[0] + z * (p[1] +
 398                  z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 399  s = q[0] + z * (q[1] + z * (q[2] + z * (q[3] + z * (q[4] + z * (q[5] + z)))));
 400   return one + z * r / s;
 401 }
 402
 403
 404 /* For x >= 8, the asymptotic expansions of qone is
 405  *      3/8 s - 105/1024 s^3 - ..., where s = 1/x.
 406  * We approximate pone by
 407  *      qone(x) = s*(0.375 + (R/S))
 408  */
 409
 410 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 411    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 412    8 <= x <= inf
 413    Peak relative error 8.3e-22 */
 414
 415 static const long double qr8[7] = {
 416   -5.691925079044209246015366919809404457380E-10L,
 417   -1.632587664706999307871963065396218379137E-7L,
 418   -1.577424682764651970003637263552027114600E-5L,
 419   -6.377627959241053914770158336842725291713E-4L,
 420   -1.087408516779972735197277149494929568768E-2L,
 421   -6.854943629378084419631926076882330494217E-2L,
 422   -1.055448290469180032312893377152490183203E-1L,
 423 };
 424 static const long double qs8[7] = {
 425   5.550982172325019811119223916998393907513E-9L,
 426   1.607188366646736068460131091130644192244E-6L,
 427   1.580792530091386496626494138334505893599E-4L,
 428   6.617859900815747303032860443855006056595E-3L,
 429   1.212840547336984859952597488863037659161E-1L,
 430   9.017885953937234900458186716154005541075E-1L,
 431   2.201114489712243262000939120146436167178E0L,
 432   /* 1.000000000000000000000000000000000000000E0L, */
 433 };
 434
 435 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 436    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 437    4.54541015625 <= x <= 8
 438    Peak relative error 4.1e-22 */
 439 static const long double qr5[7] = {
 440   -6.719134139179190546324213696633564965983E-8L,
 441   -9.467871458774950479909851595678622044140E-6L,
 442   -4.429341875348286176950914275723051452838E-4L,
 443   -8.539898021757342531563866270278505014487E-3L,
 444   -6.818691805848737010422337101409276287170E-2L,
 445   -1.964432669771684034858848142418228214855E-1L,
 446   -1.333896496989238600119596538299938520726E-1L,
 447 };
 448 static const long double qs5[7] = {
 449   6.552755584474634766937589285426911075101E-7L,
 450   9.410814032118155978663509073200494000589E-5L,
 451   4.561677087286518359461609153655021253238E-3L,
 452   9.397742096177905170800336715661091535805E-2L,
 453   8.518538116671013902180962914473967738771E-1L,
 454   3.177729183645800174212539541058292579009E0L,
 455   4.006745668510308096259753538973038902990E0L,
 456   /* 1.000000000000000000000000000000000000000E0L, */
 457 };
 458
 459 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 460    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 461    2.85711669921875 <= x <= 4.54541015625
 462    Peak relative error 2.2e-21 */
 463 static const long double qr3[7] = {
 464   -3.618746299358445926506719188614570588404E-6L,
 465   -2.951146018465419674063882650970344502798E-4L,
 466   -7.728518171262562194043409753656506795258E-3L,
 467   -8.058010968753999435006488158237984014883E-2L,
 468   -3.356232856677966691703904770937143483472E-1L,
 469   -4.858192581793118040782557808823460276452E-1L,
 470   -1.592399251246473643510898335746432479373E-1L,
 471 };
 472 static const long double qs3[7] = {
 473   3.529139957987837084554591421329876744262E-5L,
 474   2.973602667215766676998703687065066180115E-3L,
 475   8.273534546240864308494062287908662592100E-2L,
 476   9.613359842126507198241321110649974032726E-1L,
 477   4.853923697093974370118387947065402707519E0L,
 478   1.002671608961669247462020977417828796933E1L,
 479   7.028927383922483728931327850683151410267E0L,
 480   /* 1.000000000000000000000000000000000000000E0L, */
 481 };
 482
 483 /* Y1(x)cosX - J1(x)sinX = sqrt( 2/(pi x)) Q1(x),
 484    Q1(x) = z(.375 + z^2 R(z^2)), z=1/x
 485    2 <= x <= 2.85711669921875
 486    Peak relative error 6.9e-22 */
 487 static const long double qr2[7] = {
 488   -1.372751603025230017220666013816502528318E-4L,
 489   -6.879190253347766576229143006767218972834E-3L,
 490   -1.061253572090925414598304855316280077828E-1L,
 491   -6.262164224345471241219408329354943337214E-1L,
 492   -1.423149636514768476376254324731437473915E0L,
 493   -1.087955310491078933531734062917489870754E0L,
 494   -1.826821119773182847861406108689273719137E-1L,
 495 };
 496 static const long double qs2[7] = {
 497   1.338768933634451601814048220627185324007E-3L,
 498   7.071099998918497559736318523932241901810E-2L,
 499   1.200511429784048632105295629933382142221E0L,
 500   8.327301713640367079030141077172031825276E0L,
 501   2.468479301872299311658145549931764426840E1L,
 502   2.961179686096262083509383820557051621644E1L,
 503   1.201402313144305153005639494661767354977E1L,
 504  /* 1.000000000000000000000000000000000000000E0L, */
 505 };
 506
 507
 508 static long double
 509 qone (long double x)
 510 {
 511   const long double *p, *q;
 512   long double s, r, z;
 513   int32_t ix;
 514   uint32_t se, i0, i1;
 515
 516   GET_LDOUBLE_WORDS (se, i0, i1, x);
 517   ix = se & 0x7fff;
 518   /* ix >= 0x4000 for all calls to this function.  */
 519   if (ix >= 0x4002)             /* x >= 8 */
 520     {
 521       p = qr8;
 522       q = qs8;
 523     }
 524   else
 525     {
 526       i1 = (ix << 16) | (i0 >> 16);
 527       if (i1 >= 0x40019174)     /* x >= 4.54541015625 */
 528         {
 529           p = qr5;
 530           q = qs5;
 531         }
 532       else if (i1 >= 0x4000b6db)        /* x >= 2.85711669921875 */
 533         {
 534           p = qr3;
 535           q = qs3;
 536         }
 537       else      /* x >= 2 */
 538         {
 539           p = qr2;
 540           q = qs2;
 541         }
 542     }
 543   z = one / (x * x);
 544   r =
 545     p[0] + z * (p[1] +
 546                 z * (p[2] + z * (p[3] + z * (p[4] + z * (p[5] + z * p[6])))));
 547   s =
 548     q[0] + z * (q[1] +
 549                 z * (q[2] +
 550                      z * (q[3] + z * (q[4] + z * (q[5] + z * (q[6] + z))))));
 551   return (.375 + z * r / s) / x;
 552 }