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 }