00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #if defined(LIBM_SCCS) && !defined(lint)
00014 static char rcsid[] = "$NetBSD: e_j0.c,v 1.8 1995/05/10 20:45:23 jtc Exp $";
00015 #endif
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062 #include "math.h"
00063 #include "mathP.h"
00064
00065 #ifdef __STDC__
00066 static double pzero(double), qzero(double);
00067 #else
00068 static double pzero(), qzero();
00069 #endif
00070
00071 #ifdef __STDC__
00072 static const double
00073 #else
00074 static double
00075 #endif
00076 huge = 1e300,
00077 one = 1.0,
00078 invsqrtpi= 5.64189583547756279280e-01,
00079 tpi = 6.36619772367581382433e-01,
00080
00081 R02 = 1.56249999999999947958e-02,
00082 R03 = -1.89979294238854721751e-04,
00083 R04 = 1.82954049532700665670e-06,
00084 R05 = -4.61832688532103189199e-09,
00085 S01 = 1.56191029464890010492e-02,
00086 S02 = 1.16926784663337450260e-04,
00087 S03 = 5.13546550207318111446e-07,
00088 S04 = 1.16614003333790000205e-09;
00089
00090 #ifdef __STDC__
00091 static const double zero = 0.0;
00092 #else
00093 static double zero = 0.0;
00094 #endif
00095
00096 #ifdef __STDC__
00097 double __ieee754_j0(double x)
00098 #else
00099 double __ieee754_j0(x)
00100 double x;
00101 #endif
00102 {
00103 double z, s,c,ss,cc,r,u,v;
00104 int32_t hx,ix;
00105
00106 GET_HIGH_WORD(hx,x);
00107 ix = hx&0x7fffffff;
00108 if(ix>=0x7ff00000) return one/(x*x);
00109 x = fabs(x);
00110 if(ix >= 0x40000000) {
00111 s = sin(x);
00112 c = cos(x);
00113 ss = s-c;
00114 cc = s+c;
00115 if(ix<0x7fe00000) {
00116 z = -cos(x+x);
00117 if ((s*c)<zero) cc = z/ss;
00118 else ss = z/cc;
00119 }
00120
00121
00122
00123
00124 if(ix>0x48000000) z = (invsqrtpi*cc)/sqrt(x);
00125 else {
00126 u = pzero(x); v = qzero(x);
00127 z = invsqrtpi*(u*cc-v*ss)/sqrt(x);
00128 }
00129 return z;
00130 }
00131 if(ix<0x3f200000) {
00132 if(huge+x>one) {
00133 if(ix<0x3e400000) return one;
00134 else return one - 0.25*x*x;
00135 }
00136 }
00137 z = x*x;
00138 r = z*(R02+z*(R03+z*(R04+z*R05)));
00139 s = one+z*(S01+z*(S02+z*(S03+z*S04)));
00140 if(ix < 0x3FF00000) {
00141 return one + z*(-0.25+(r/s));
00142 } else {
00143 u = 0.5*x;
00144 return((one+u)*(one-u)+z*(r/s));
00145 }
00146 }
00147
00148 #ifdef __STDC__
00149 static const double
00150 #else
00151 static double
00152 #endif
00153 u00 = -7.38042951086872317523e-02,
00154 u01 = 1.76666452509181115538e-01,
00155 u02 = -1.38185671945596898896e-02,
00156 u03 = 3.47453432093683650238e-04,
00157 u04 = -3.81407053724364161125e-06,
00158 u05 = 1.95590137035022920206e-08,
00159 u06 = -3.98205194132103398453e-11,
00160 v01 = 1.27304834834123699328e-02,
00161 v02 = 7.60068627350353253702e-05,
00162 v03 = 2.59150851840457805467e-07,
00163 v04 = 4.41110311332675467403e-10;
00164
00165 #ifdef __STDC__
00166 double __ieee754_y0(double x)
00167 #else
00168 double __ieee754_y0(x)
00169 double x;
00170 #endif
00171 {
00172 double z, s,c,ss,cc,u,v;
00173 int32_t hx,ix,lx;
00174
00175 EXTRACT_WORDS(hx,lx,x);
00176 ix = 0x7fffffff&hx;
00177
00178 if(ix>=0x7ff00000) return one/(x+x*x);
00179 if((ix|lx)==0) return -one/zero;
00180 if(hx<0) return zero/zero;
00181 if(ix >= 0x40000000) {
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193 s = sin(x);
00194 c = cos(x);
00195 ss = s-c;
00196 cc = s+c;
00197
00198
00199
00200
00201 if(ix<0x7fe00000) {
00202 z = -cos(x+x);
00203 if ((s*c)<zero) cc = z/ss;
00204 else ss = z/cc;
00205 }
00206 if(ix>0x48000000) z = (invsqrtpi*ss)/sqrt(x);
00207 else {
00208 u = pzero(x); v = qzero(x);
00209 z = invsqrtpi*(u*ss+v*cc)/sqrt(x);
00210 }
00211 return z;
00212 }
00213 if(ix<=0x3e400000) {
00214 return(u00 + tpi*__ieee754_log(x));
00215 }
00216 z = x*x;
00217 u = u00+z*(u01+z*(u02+z*(u03+z*(u04+z*(u05+z*u06)))));
00218 v = one+z*(v01+z*(v02+z*(v03+z*v04)));
00219 return(u/v + tpi*(__ieee754_j0(x)*__ieee754_log(x)));
00220 }
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231 #ifdef __STDC__
00232 static const double pR8[6] = {
00233 #else
00234 static double pR8[6] = {
00235 #endif
00236 0.00000000000000000000e+00,
00237 -7.03124999999900357484e-02,
00238 -8.08167041275349795626e+00,
00239 -2.57063105679704847262e+02,
00240 -2.48521641009428822144e+03,
00241 -5.25304380490729545272e+03,
00242 };
00243 #ifdef __STDC__
00244 static const double pS8[5] = {
00245 #else
00246 static double pS8[5] = {
00247 #endif
00248 1.16534364619668181717e+02,
00249 3.83374475364121826715e+03,
00250 4.05978572648472545552e+04,
00251 1.16752972564375915681e+05,
00252 4.76277284146730962675e+04,
00253 };
00254
00255 #ifdef __STDC__
00256 static const double pR5[6] = {
00257 #else
00258 static double pR5[6] = {
00259 #endif
00260 -1.14125464691894502584e-11,
00261 -7.03124940873599280078e-02,
00262 -4.15961064470587782438e+00,
00263 -6.76747652265167261021e+01,
00264 -3.31231299649172967747e+02,
00265 -3.46433388365604912451e+02,
00266 };
00267 #ifdef __STDC__
00268 static const double pS5[5] = {
00269 #else
00270 static double pS5[5] = {
00271 #endif
00272 6.07539382692300335975e+01,
00273 1.05125230595704579173e+03,
00274 5.97897094333855784498e+03,
00275 9.62544514357774460223e+03,
00276 2.40605815922939109441e+03,
00277 };
00278
00279 #ifdef __STDC__
00280 static const double pR3[6] = {
00281 #else
00282 static double pR3[6] = {
00283 #endif
00284 -2.54704601771951915620e-09,
00285 -7.03119616381481654654e-02,
00286 -2.40903221549529611423e+00,
00287 -2.19659774734883086467e+01,
00288 -5.80791704701737572236e+01,
00289 -3.14479470594888503854e+01,
00290 };
00291 #ifdef __STDC__
00292 static const double pS3[5] = {
00293 #else
00294 static double pS3[5] = {
00295 #endif
00296 3.58560338055209726349e+01,
00297 3.61513983050303863820e+02,
00298 1.19360783792111533330e+03,
00299 1.12799679856907414432e+03,
00300 1.73580930813335754692e+02,
00301 };
00302
00303 #ifdef __STDC__
00304 static const double pR2[6] = {
00305 #else
00306 static double pR2[6] = {
00307 #endif
00308 -8.87534333032526411254e-08,
00309 -7.03030995483624743247e-02,
00310 -1.45073846780952986357e+00,
00311 -7.63569613823527770791e+00,
00312 -1.11931668860356747786e+01,
00313 -3.23364579351335335033e+00,
00314 };
00315 #ifdef __STDC__
00316 static const double pS2[5] = {
00317 #else
00318 static double pS2[5] = {
00319 #endif
00320 2.22202997532088808441e+01,
00321 1.36206794218215208048e+02,
00322 2.70470278658083486789e+02,
00323 1.53875394208320329881e+02,
00324 1.46576176948256193810e+01,
00325 };
00326
00327 #ifdef __STDC__
00328 static double pzero(double x)
00329 #else
00330 static double pzero(x)
00331 double x;
00332 #endif
00333 {
00334 #ifdef __STDC__
00335 const double *p = 0,*q = 0;
00336 #else
00337 double *p = 0,*q = 0;
00338 #endif
00339 double z,r,s;
00340 int32_t ix;
00341 GET_HIGH_WORD(ix,x);
00342 ix &= 0x7fffffff;
00343 if(ix>=0x40200000) {p = pR8; q= pS8;}
00344 else if(ix>=0x40122E8B){p = pR5; q= pS5;}
00345 else if(ix>=0x4006DB6D){p = pR3; q= pS3;}
00346 else if(ix>=0x40000000){p = pR2; q= pS2;}
00347 z = one/(x*x);
00348 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
00349 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
00350 return one+ r/s;
00351 }
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363 #ifdef __STDC__
00364 static const double qR8[6] = {
00365 #else
00366 static double qR8[6] = {
00367 #endif
00368 0.00000000000000000000e+00,
00369 7.32421874999935051953e-02,
00370 1.17682064682252693899e+01,
00371 5.57673380256401856059e+02,
00372 8.85919720756468632317e+03,
00373 3.70146267776887834771e+04,
00374 };
00375 #ifdef __STDC__
00376 static const double qS8[6] = {
00377 #else
00378 static double qS8[6] = {
00379 #endif
00380 1.63776026895689824414e+02,
00381 8.09834494656449805916e+03,
00382 1.42538291419120476348e+05,
00383 8.03309257119514397345e+05,
00384 8.40501579819060512818e+05,
00385 -3.43899293537866615225e+05,
00386 };
00387
00388 #ifdef __STDC__
00389 static const double qR5[6] = {
00390 #else
00391 static double qR5[6] = {
00392 #endif
00393 1.84085963594515531381e-11,
00394 7.32421766612684765896e-02,
00395 5.83563508962056953777e+00,
00396 1.35111577286449829671e+02,
00397 1.02724376596164097464e+03,
00398 1.98997785864605384631e+03,
00399 };
00400 #ifdef __STDC__
00401 static const double qS5[6] = {
00402 #else
00403 static double qS5[6] = {
00404 #endif
00405 8.27766102236537761883e+01,
00406 2.07781416421392987104e+03,
00407 1.88472887785718085070e+04,
00408 5.67511122894947329769e+04,
00409 3.59767538425114471465e+04,
00410 -5.35434275601944773371e+03,
00411 };
00412
00413 #ifdef __STDC__
00414 static const double qR3[6] = {
00415 #else
00416 static double qR3[6] = {
00417 #endif
00418 4.37741014089738620906e-09,
00419 7.32411180042911447163e-02,
00420 3.34423137516170720929e+00,
00421 4.26218440745412650017e+01,
00422 1.70808091340565596283e+02,
00423 1.66733948696651168575e+02,
00424 };
00425 #ifdef __STDC__
00426 static const double qS3[6] = {
00427 #else
00428 static double qS3[6] = {
00429 #endif
00430 4.87588729724587182091e+01,
00431 7.09689221056606015736e+02,
00432 3.70414822620111362994e+03,
00433 6.46042516752568917582e+03,
00434 2.51633368920368957333e+03,
00435 -1.49247451836156386662e+02,
00436 };
00437
00438 #ifdef __STDC__
00439 static const double qR2[6] = {
00440 #else
00441 static double qR2[6] = {
00442 #endif
00443 1.50444444886983272379e-07,
00444 7.32234265963079278272e-02,
00445 1.99819174093815998816e+00,
00446 1.44956029347885735348e+01,
00447 3.16662317504781540833e+01,
00448 1.62527075710929267416e+01,
00449 };
00450 #ifdef __STDC__
00451 static const double qS2[6] = {
00452 #else
00453 static double qS2[6] = {
00454 #endif
00455 3.03655848355219184498e+01,
00456 2.69348118608049844624e+02,
00457 8.44783757595320139444e+02,
00458 8.82935845112488550512e+02,
00459 2.12666388511798828631e+02,
00460 -5.31095493882666946917e+00,
00461 };
00462
00463 #ifdef __STDC__
00464 static double qzero(double x)
00465 #else
00466 static double qzero(x)
00467 double x;
00468 #endif
00469 {
00470 #ifdef __STDC__
00471 const double *p = 0,*q = 0;
00472 #else
00473 double *p = 0,*q = 0;
00474 #endif
00475 double s,r,z;
00476 int32_t ix;
00477 GET_HIGH_WORD(ix,x);
00478 ix &= 0x7fffffff;
00479 if(ix>=0x40200000) {p = qR8; q= qS8;}
00480 else if(ix>=0x40122E8B){p = qR5; q= qS5;}
00481 else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
00482 else if(ix>=0x40000000){p = qR2; q= qS2;}
00483 z = one/(x*x);
00484 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
00485 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
00486 return (-.125 + r/s)/x;
00487 }