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 __j0 = j0
  40 #pragma weak __y0 = y0
  41 
  42 #include "libm.h"
  43 #include "libm_protos.h"
  44 #include <math.h>
  45 #include <values.h>
  46 
  47 #define GENERIC double
  48 static const GENERIC
  49 zero    = 0.0,
  50 small   = 1.0e-5,
  51 tiny    = 1.0e-18,
  52 one     = 1.0,
  53 eight   = 8.0,
  54 invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
  55 tpi     = 0.636619772367581343075535053490057448;
  56 
  57 static GENERIC pzero(GENERIC), qzero(GENERIC);










  58 static const GENERIC r0[4] = {  /* [1.e-5, 1.28] */
  59         -2.500000000000003622131880894830476755537e-0001,
  60         1.095597547334830263234433855932375353303e-0002,
  61         -1.819734750463320921799187258987098087697e-0004,
  62         9.977001946806131657544212501069893930846e-0007,
  63 };

  64 static const GENERIC s0[4] = {  /* [1.e-5, 1.28] */
  65         1.0,
  66         1.867609810662950169966782360588199673741e-0002,
  67         1.590389206181565490878430827706972074208e-0004,
  68         6.520867386742583632375520147714499522721e-0007,
  69 };

  70 static const GENERIC r1[9] = {  /* [1.28,8] */
  71         9.999999999999999942156495584397047660949e-0001,
  72         -2.389887722731319130476839836908143731281e-0001,
  73         1.293359476138939027791270393439493640570e-0002,
  74         -2.770985642343140122168852400228563364082e-0004,
  75         2.905241575772067678086738389169625218912e-0006,
  76         -1.636846356264052597969042009265043251279e-0008,
  77         5.072306160724884775085431059052611737827e-0011,
  78         -8.187060730684066824228914775146536139112e-0014,
  79         5.422219326959949863954297860723723423842e-0017,
  80 };

  81 static const GENERIC s1[9] = {  /* [1.28,8] */
  82         1.0,
  83         1.101122772686807702762104741932076228349e-0002,
  84         6.140169310641649223411427764669143978228e-0005,
  85         2.292035877515152097976946119293215705250e-0007,
  86         6.356910426504644334558832036362219583789e-0010,
  87         1.366626326900219555045096999553948891401e-0012,
  88         2.280399586866739522891837985560481180088e-0015,
  89         2.801559820648939665270492520004836611187e-0018,
  90         2.073101088320349159764410261466350732968e-0021,
  91 };
  92 
  93 GENERIC
  94 j0(GENERIC x) {

  95         GENERIC z, s, c, ss, cc, r, u, v, ox;
  96         int i;
  97 
  98         if (isnan(x))
  99                 return (x*x);   /* + -> * for Cheetah */

 100         ox = x;
 101         x = fabs(x);

 102         if (x > 8.0) {
 103                 if (!finite(x))
 104                         return (zero);

 105                 s = sin(x);
 106                 c = cos(x);


 107         /*
 108          * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
 109          * where x0 = x-pi/4
 110          *      Better formula:
 111          *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
 112          *                      = 1/sqrt(2) * (cos(x) + sin(x))
 113          *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
 114          *                      = 1/sqrt(2) * (sin(x) - cos(x))
 115          * To avoid cancellation, use
 116          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 117          * to compute the worse one.
 118          */

 119                 if (x > 8.9e307) {   /* x+x may overflow */
 120                         ss = s-c;
 121                         cc = s+c;
 122                 } else if (signbit(s) != signbit(c)) {
 123                         ss = s - c;
 124                         cc = -cos(x+x)/ss;
 125                 } else {
 126                         cc = s + c;
 127                         ss = -cos(x+x)/cc;
 128                 }

 129         /*
 130          * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 131          * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 132          */
 133                 if (x > 1.0e40) z = (invsqrtpi*cc)/sqrt(x);
 134                 else {
 135                     u = pzero(x); v = qzero(x);
 136                     z = invsqrtpi*(u*cc-v*ss)/sqrt(x);


 137                 }

 138         /* force to pass SVR4 even the result is wrong (sign) */
 139                 if (x > X_TLOSS)
 140                     return (_SVID_libm_err(ox, z, 34));
 141                 else
 142                     return (z);
 143         }

 144         if (x <= small) {
 145             if (x <= tiny)
 146                         return (one-x);
 147             else
 148                         return (one-x*x*0.25);
 149         }
 150         z = x*x;


 151         if (x <= 1.28) {
 152             r =  r0[0]+z*(r0[1]+z*(r0[2]+z*r0[3]));
 153             s =  s0[0]+z*(s0[1]+z*(s0[2]+z*s0[3]));
 154             return (one + z*(r/s));
 155         } else {
 156             for (r = r1[8], s = s1[8], i = 7; i >= 0; i--) {
 157                 r = r*z + r1[i];
 158                 s = s*z + s1[i];
 159             }
 160             return (r/s);

 161         }
 162 }
 163 
 164 static const GENERIC u0[13] = {
 165         -7.380429510868722526754723020704317641941e-0002,
 166         1.772607102684869924301459663049874294814e-0001,
 167         -1.524370666542713828604078090970799356306e-0002,
 168         4.650819100693891757143771557629924591915e-0004,
 169         -7.125768872339528975036316108718239946022e-0006,
 170         6.411017001656104598327565004771515257146e-0008,
 171         -3.694275157433032553021246812379258781665e-0010,
 172         1.434364544206266624252820889648445263842e-0012,
 173         -3.852064731859936455895036286874139896861e-0015,
 174         7.182052899726138381739945881914874579696e-0018,
 175         -9.060556574619677567323741194079797987200e-0021,
 176         7.124435467408860515265552217131230511455e-0024,
 177         -2.709726774636397615328813121715432044771e-0027,
 178 };

 179 static const GENERIC v0[5] = {
 180         1.0,
 181         4.678678931512549002587702477349214886475e-0003,
 182         9.486828955529948534822800829497565178985e-0006,
 183         1.001495929158861646659010844136682454906e-0008,
 184         4.725338116256021660204443235685358593611e-0012,
 185 };
 186 
 187 GENERIC
 188 y0(GENERIC x) {

 189         GENERIC z, /* d, */ s, c, ss, cc, u, v;
 190         int i;
 191 
 192         if (isnan(x))
 193                 return (x*x);   /* + -> * for Cheetah */

 194         if (x <= zero) {
 195                 if (x == zero)
 196                     /* d= -one/(x-x); */
 197                     return (_SVID_libm_err(x, x, 8));
 198                 else
 199                     /* d = zero/(x-x); */
 200                     return (_SVID_libm_err(x, x, 9));
 201         }

 202         if (x > 8.0) {
 203                 if (!finite(x))
 204                         return (zero);

 205                 s = sin(x);
 206                 c = cos(x);


 207         /*
 208          * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
 209          * where x0 = x-pi/4
 210          *      Better formula:
 211          *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
 212          *                      = 1/sqrt(2) * (cos(x) + sin(x))
 213          *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
 214          *                      = 1/sqrt(2) * (sin(x) - cos(x))
 215          * To avoid cancellation, use
 216          *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 217          * to compute the worse one.
 218          */

 219                 if (x > 8.9e307) {   /* x+x may overflow */
 220                         ss = s-c;
 221                         cc = s+c;
 222                 } else if (signbit(s) != signbit(c)) {
 223                         ss = s - c;
 224                         cc = -cos(x+x)/ss;
 225                 } else {
 226                         cc = s + c;
 227                         ss = -cos(x+x)/cc;
 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.0e40)
 234                     z = (invsqrtpi*ss)/sqrt(x);
 235                 else
 236                     z =  invsqrtpi*(pzero(x)*ss+qzero(x)*cc)/sqrt(x);


 237                 if (x > X_TLOSS)
 238                     return (_SVID_libm_err(x, z, 35));
 239                 else
 240                     return (z);
 241 
 242         }
 243         if (x <= tiny) {
 244             return (u0[0] + tpi*log(x));
 245         }
 246         z = x*x;
 247         for (u = u0[12], i = 11; i >= 0; i--) u = u*z + u0[i];
 248         v = v0[0]+z*(v0[1]+z*(v0[2]+z*(v0[3]+z*v0[4])));
 249         return (u/v + tpi*(j0(x)*log(x)));







 250 }
 251 
 252 static const GENERIC pr[7] = {  /* [8 -- inf]  pzero 6550 */
 253         .4861344183386052721391238447e5,
 254         .1377662549407112278133438945e6,
 255         .1222466364088289731869114004e6,
 256         .4107070084315176135583353374e5,
 257         .5026073801860637125889039915e4,
 258         .1783193659125479654541542419e3,
 259         .88010344055383421691677564e0,
 260 };

 261 static const GENERIC ps[7] = {  /* [8 -- inf] pzero 6550 */
 262         .4861344183386052721414037058e5,
 263         .1378196632630384670477582699e6,
 264         .1223967185341006542748936787e6,
 265         .4120150243795353639995862617e5,
 266         .5068271181053546392490184353e4,
 267         .1829817905472769960535671664e3,
 268         1.0,
 269 };
 270 static const GENERIC huge    = 1.0e10;
 271 

 272 static GENERIC
 273 pzero(GENERIC x) {

 274         GENERIC s, r, t, z;
 275         int i;

 276         if (x > huge)
 277                 return (one);
 278         t = eight/x; z = t*t;
 279         r = pr[5]+z*pr[6];
 280         s = ps[5]+z;



 281         for (i = 4; i >= 0; i--) {
 282             r = r*z + pr[i];
 283             s = s*z + ps[i];
 284         }
 285         return (r/s);

 286 }
 287 
 288 static const GENERIC qr[7] = {  /* [8 -- inf]  qzero 6950 */
 289         -.1731210995701068539185611951e3,
 290         -.5522559165936166961235240613e3,
 291         -.5604935606637346590614529613e3,
 292         -.2200430300226009379477365011e3,
 293         -.323869355375648849771296746e2,
 294         -.14294979207907956223499258e1,
 295         -.834690374102384988158918e-2,
 296 };

 297 static const GENERIC qs[7] = {  /* [8 -- inf] qzero 6950 */
 298         .1107975037248683865326709645e5,
 299         .3544581680627082674651471873e5,
 300         .3619118937918394132179019059e5,
 301         .1439895563565398007471485822e5,
 302         .2190277023344363955930226234e4,
 303         .106695157020407986137501682e3,
 304         1.0,
 305 };
 306 
 307 static GENERIC
 308 qzero(GENERIC x) {

 309         GENERIC s, r, t, z;
 310         int i;

 311         if (x > huge)
 312                 return (-0.125/x);
 313         t = eight/x; z = t*t;
 314         r = qr[5]+z*qr[6];
 315         s = qs[5]+z;



 316         for (i = 4; i >= 0; i--) {
 317             r = r*z + qr[i];
 318             s = s*z + qs[i];
 319         }
 320         return (t*(r/s));

 321 }


   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 __j0 = j0
  41 #pragma weak __y0 = y0
  42 
  43 #include "libm.h"
  44 #include "libm_protos.h"
  45 #include <math.h>
  46 #include <values.h>
  47 
  48 #define GENERIC double








  49 
  50 static const GENERIC zero = 0.0,
  51         small = 1.0e-5,
  52         tiny = 1.0e-18,
  53         one = 1.0,
  54         eight = 8.0,
  55         invsqrtpi = 5.641895835477562869480794515607725858441e-0001,
  56         tpi = 0.636619772367581343075535053490057448;
  57 
  58 static GENERIC pzero(GENERIC);
  59 static GENERIC qzero(GENERIC);
  60 
  61 static const GENERIC r0[4] = {          /* [1.e-5, 1.28] */
  62         -2.500000000000003622131880894830476755537e-0001,
  63         1.095597547334830263234433855932375353303e-0002,
  64         -1.819734750463320921799187258987098087697e-0004,
  65         9.977001946806131657544212501069893930846e-0007,
  66 };
  67 
  68 static const GENERIC s0[4] = {          /* [1.e-5, 1.28] */
  69         1.0,
  70         1.867609810662950169966782360588199673741e-0002,
  71         1.590389206181565490878430827706972074208e-0004,
  72         6.520867386742583632375520147714499522721e-0007,
  73 };
  74 
  75 static const GENERIC r1[9] = {          /* [1.28,8] */
  76         9.999999999999999942156495584397047660949e-0001,
  77         -2.389887722731319130476839836908143731281e-0001,
  78         1.293359476138939027791270393439493640570e-0002,
  79         -2.770985642343140122168852400228563364082e-0004,
  80         2.905241575772067678086738389169625218912e-0006,
  81         -1.636846356264052597969042009265043251279e-0008,
  82         5.072306160724884775085431059052611737827e-0011,
  83         -8.187060730684066824228914775146536139112e-0014,
  84         5.422219326959949863954297860723723423842e-0017,
  85 };
  86 
  87 static const GENERIC s1[9] = {          /* [1.28,8] */
  88         1.0,
  89         1.101122772686807702762104741932076228349e-0002,
  90         6.140169310641649223411427764669143978228e-0005,
  91         2.292035877515152097976946119293215705250e-0007,
  92         6.356910426504644334558832036362219583789e-0010,
  93         1.366626326900219555045096999553948891401e-0012,
  94         2.280399586866739522891837985560481180088e-0015,
  95         2.801559820648939665270492520004836611187e-0018,
  96         2.073101088320349159764410261466350732968e-0021,
  97 };
  98 
  99 GENERIC
 100 j0(GENERIC x)
 101 {
 102         GENERIC z, s, c, ss, cc, r, u, v, ox;
 103         int i;
 104 
 105         if (isnan(x))
 106                 return (x * x);         /* + -> * for Cheetah */
 107 
 108         ox = x;
 109         x = fabs(x);
 110 
 111         if (x > 8.0) {
 112                 if (!finite(x))
 113                         return (zero);
 114 
 115                 s = sin(x);
 116                 c = cos(x);
 117 
 118                 /* BEGIN CSTYLED */
 119                 /*
 120                  * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
 121                  * where x0 = x-pi/4
 122                  *      Better formula:
 123                  *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
 124                  *                      = 1/sqrt(2) * (cos(x) + sin(x))
 125                  *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
 126                  *                      = 1/sqrt(2) * (sin(x) - cos(x))
 127                  * To avoid cancellation, use
 128                  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 129                  * to compute the worse one.
 130                  */
 131                 /* END CSTYLED */
 132                 if (x > 8.9e307) {   /* x+x may overflow */
 133                         ss = s - c;
 134                         cc = s + c;
 135                 } else if (signbit(s) != signbit(c)) {
 136                         ss = s - c;
 137                         cc = -cos(x + x) / ss;
 138                 } else {
 139                         cc = s + c;
 140                         ss = -cos(x + x) / cc;
 141                 }
 142 
 143                 /*
 144                  * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
 145                  * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
 146                  */
 147                 if (x > 1.0e40) {
 148                         z = (invsqrtpi * cc) / sqrt(x);
 149                 } else {
 150                         u = pzero(x);
 151                         v = qzero(x);
 152                         z = invsqrtpi * (u * cc - v * ss) / sqrt(x);
 153                 }
 154 
 155                 /* force to pass SVR4 even the result is wrong (sign) */
 156                 if (x > X_TLOSS)
 157                         return (_SVID_libm_err(ox, z, 34));
 158                 else
 159                         return (z);
 160         }
 161 
 162         if (x <= small) {
 163                 if (x <= tiny)
 164                         return (one - x);
 165                 else
 166                         return (one - x * x * 0.25);
 167         }
 168 
 169         z = x * x;
 170 
 171         if (x <= 1.28) {
 172                 r = r0[0] + z * (r0[1] + z * (r0[2] + z * r0[3]));
 173                 s = s0[0] + z * (s0[1] + z * (s0[2] + z * s0[3]));
 174                 return (one + z * (r / s));
 175         } else {
 176                 for (r = r1[8], s = s1[8], i = 7; i >= 0; i--) {
 177                         r = r * z + r1[i];
 178                         s = s * z + s1[i];
 179                 }
 180 
 181                 return (r / s);
 182         }
 183 }
 184 
 185 static const GENERIC u0[13] = {
 186         -7.380429510868722526754723020704317641941e-0002,
 187         1.772607102684869924301459663049874294814e-0001,
 188         -1.524370666542713828604078090970799356306e-0002,
 189         4.650819100693891757143771557629924591915e-0004,
 190         -7.125768872339528975036316108718239946022e-0006,
 191         6.411017001656104598327565004771515257146e-0008,
 192         -3.694275157433032553021246812379258781665e-0010,
 193         1.434364544206266624252820889648445263842e-0012,
 194         -3.852064731859936455895036286874139896861e-0015,
 195         7.182052899726138381739945881914874579696e-0018,
 196         -9.060556574619677567323741194079797987200e-0021,
 197         7.124435467408860515265552217131230511455e-0024,
 198         -2.709726774636397615328813121715432044771e-0027,
 199 };
 200 
 201 static const GENERIC v0[5] = {
 202         1.0,
 203         4.678678931512549002587702477349214886475e-0003,
 204         9.486828955529948534822800829497565178985e-0006,
 205         1.001495929158861646659010844136682454906e-0008,
 206         4.725338116256021660204443235685358593611e-0012,
 207 };
 208 
 209 GENERIC
 210 y0(GENERIC x)
 211 {
 212         GENERIC z, /* d, */ s, c, ss, cc, u, v;
 213         int i;
 214 
 215         if (isnan(x))
 216                 return (x * x);         /* + -> * for Cheetah */
 217 
 218         if (x <= zero) {
 219                 if (x == zero)
 220                         /* d= -one/(x-x); */
 221                         return (_SVID_libm_err(x, x, 8));
 222                 else
 223                         /* d = zero/(x-x); */
 224                         return (_SVID_libm_err(x, x, 9));
 225         }
 226 
 227         if (x > 8.0) {
 228                 if (!finite(x))
 229                         return (zero);
 230 
 231                 s = sin(x);
 232                 c = cos(x);
 233 
 234                 /* BEGIN CSTYLED */
 235                 /*
 236                  * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
 237                  * where x0 = x-pi/4
 238                  *      Better formula:
 239                  *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
 240                  *                      = 1/sqrt(2) * (cos(x) + sin(x))
 241                  *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
 242                  *                      = 1/sqrt(2) * (sin(x) - cos(x))
 243                  * To avoid cancellation, use
 244                  *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
 245                  * to compute the worse one.
 246                  */
 247                 /* END CSTYLED */
 248                 if (x > 8.9e307) {   /* x+x may overflow */
 249                         ss = s - c;
 250                         cc = s + c;
 251                 } else if (signbit(s) != signbit(c)) {
 252                         ss = s - c;
 253                         cc = -cos(x + x) / ss;
 254                 } else {
 255                         cc = s + c;
 256                         ss = -cos(x + x) / cc;
 257                 }
 258 
 259                 /*
 260                  * j0(x) = 1/sqrt(pi*x) * (P(0,x)*cc - Q(0,x)*ss)
 261                  * y0(x) = 1/sqrt(pi*x) * (P(0,x)*ss + Q(0,x)*cc)
 262                  */
 263                 if (x > 1.0e40)
 264                         z = (invsqrtpi * ss) / sqrt(x);
 265                 else
 266                         z = invsqrtpi * (pzero(x) * ss + qzero(x) * cc) /
 267                             sqrt(x);
 268 
 269                 if (x > X_TLOSS)
 270                         return (_SVID_libm_err(x, z, 35));
 271                 else
 272                         return (z);




 273         }
 274 
 275         if (x <= tiny)
 276                 return (u0[0] + tpi * log(x));
 277 
 278         z = x * x;
 279 
 280         for (u = u0[12], i = 11; i >= 0; i--)
 281                 u = u * z + u0[i];
 282 
 283         v = v0[0] + z * (v0[1] + z * (v0[2] + z * (v0[3] + z * v0[4])));
 284         return (u / v + tpi * (j0(x) * log(x)));
 285 }
 286 
 287 static const GENERIC pr[7] = {          /* [8 -- inf]  pzero 6550 */
 288         .4861344183386052721391238447e5, .1377662549407112278133438945e6,
 289         .1222466364088289731869114004e6, .4107070084315176135583353374e5,
 290         .5026073801860637125889039915e4, .1783193659125479654541542419e3,



 291         .88010344055383421691677564e0,
 292 };
 293 
 294 static const GENERIC ps[7] = {          /* [8 -- inf] pzero 6550 */
 295         .4861344183386052721414037058e5, .1378196632630384670477582699e6,
 296         .1223967185341006542748936787e6, .4120150243795353639995862617e5,
 297         .5068271181053546392490184353e4, .1829817905472769960535671664e3,



 298         1.0,
 299 };

 300 
 301 static const GENERIC huge = 1.0e10;
 302 static GENERIC
 303 pzero(GENERIC x)
 304 {
 305         GENERIC s, r, t, z;
 306         int i;
 307 
 308         if (x > huge)
 309                 return (one);
 310 
 311         t = eight / x;
 312         z = t * t;
 313         r = pr[5] + z * pr[6];
 314         s = ps[5] + z;
 315 
 316         for (i = 4; i >= 0; i--) {
 317                 r = r * z + pr[i];
 318                 s = s * z + ps[i];
 319         }
 320 
 321         return (r / s);
 322 }
 323 
 324 static const GENERIC qr[7] = {          /* [8 -- inf]  qzero 6950 */
 325         -.1731210995701068539185611951e3, -.5522559165936166961235240613e3,
 326         -.5604935606637346590614529613e3, -.2200430300226009379477365011e3,
 327         -.323869355375648849771296746e2, -.14294979207907956223499258e1,



 328         -.834690374102384988158918e-2,
 329 };
 330 
 331 static const GENERIC qs[7] = {          /* [8 -- inf] qzero 6950 */
 332         .1107975037248683865326709645e5, .3544581680627082674651471873e5,
 333         .3619118937918394132179019059e5, .1439895563565398007471485822e5,
 334         .2190277023344363955930226234e4, .106695157020407986137501682e3,



 335         1.0,
 336 };
 337 
 338 static GENERIC
 339 qzero(GENERIC x)
 340 {
 341         GENERIC s, r, t, z;
 342         int i;
 343 
 344         if (x > huge)
 345                 return (-0.125 / x);
 346 
 347         t = eight / x;
 348         z = t * t;
 349         r = qr[5] + z * qr[6];
 350         s = qs[5] + z;
 351 
 352         for (i = 4; i >= 0; i--) {
 353                 r = r * z + qr[i];
 354                 s = s * z + qs[i];
 355         }
 356 
 357         return (t * (r / s));
 358 }