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_j1.c,v 1.8 1995/05/10 20:45:27 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 pone(
double),
qone(
double);
00067
#else
00068
static double pone(),
qone();
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 r00 = -6.25000000000000000000e-02,
00082 r01 = 1.40705666955189706048e-03,
00083 r02 = -1.59955631084035597520e-05,
00084 r03 = 4.96727999609584448412e-08,
00085 s01 = 1.91537599538363460805e-02,
00086 s02 = 1.85946785588630915560e-04,
00087 s03 = 1.17718464042623683263e-06,
00088 s04 = 5.04636257076217042715e-09,
00089 s05 = 1.23542274426137913908e-11;
00090
00091
#ifdef __STDC__
00092
static const double zero = 0.0;
00093
#else
00094 static double zero = 0.0;
00095
#endif
00096
00097
#ifdef __STDC__
00098
double __ieee754_j1(
double x)
00099 #
else
00100 double __ieee754_j1(x)
00101 double x;
00102 #endif
00103 {
00104
double z, s,c,ss,cc,r,u,v,y;
00105 int32_t hx,ix;
00106
00107
GET_HIGH_WORD(hx,x);
00108 ix = hx&0x7fffffff;
00109
if(ix>=0x7ff00000)
return one/x;
00110 y =
fabs(x);
00111
if(ix >= 0x40000000) {
00112 s =
sin(y);
00113 c =
cos(y);
00114 ss = -s-c;
00115 cc = s-c;
00116
if(ix<0x7fe00000) {
00117 z =
cos(y+y);
00118
if ((s*c)>
zero) cc = z/ss;
00119
else ss = z/cc;
00120 }
00121
00122
00123
00124
00125
if(ix>0x48000000) z = (invsqrtpi*cc)/
sqrt(y);
00126
else {
00127 u =
pone(y); v =
qone(y);
00128 z = invsqrtpi*(u*cc-v*ss)/
sqrt(y);
00129 }
00130
if(hx<0)
return -z;
00131
else return z;
00132 }
00133
if(ix<0x3e400000) {
00134
if(
huge+x>one)
return 0.5*x;
00135 }
00136 z = x*x;
00137 r = z*(r00+z*(r01+z*(r02+z*r03)));
00138 s = one+z*(s01+z*(s02+z*(s03+z*(s04+z*s05))));
00139 r *= x;
00140
return(x*0.5+r/s);
00141 }
00142
00143
#ifdef __STDC__
00144
static const double U0[5] = {
00145
#else
00146 static double U0[5] = {
00147
#endif
00148
-1.96057090646238940668e-01,
00149 5.04438716639811282616e-02,
00150 -1.91256895875763547298e-03,
00151 2.35252600561610495928e-05,
00152 -9.19099158039878874504e-08,
00153 };
00154
#ifdef __STDC__
00155
static const double V0[5] = {
00156
#else
00157 static double V0[5] = {
00158
#endif
00159
1.99167318236649903973e-02,
00160 2.02552581025135171496e-04,
00161 1.35608801097516229404e-06,
00162 6.22741452364621501295e-09,
00163 1.66559246207992079114e-11,
00164 };
00165
00166
#ifdef __STDC__
00167
double __ieee754_y1(
double x)
00168 #
else
00169 double __ieee754_y1(x)
00170 double x;
00171 #endif
00172 {
00173
double z, s,c,ss,cc,u,v;
00174 int32_t hx,ix,lx;
00175
00176
EXTRACT_WORDS(hx,lx,x);
00177 ix = 0x7fffffff&hx;
00178
00179
if(ix>=0x7ff00000)
return one/(x+x*x);
00180
if((ix|lx)==0)
return -one/
zero;
00181
if(hx<0)
return zero/
zero;
00182
if(ix >= 0x40000000) {
00183 s =
sin(x);
00184 c =
cos(x);
00185 ss = -s-c;
00186 cc = s-c;
00187
if(ix<0x7fe00000) {
00188 z =
cos(x+x);
00189
if ((s*c)>
zero) cc = z/ss;
00190
else ss = z/cc;
00191 }
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
if(ix>0x48000000) z = (invsqrtpi*ss)/
sqrt(x);
00204
else {
00205 u =
pone(x); v =
qone(x);
00206 z = invsqrtpi*(u*ss+v*cc)/
sqrt(x);
00207 }
00208
return z;
00209 }
00210
if(ix<=0x3c900000) {
00211
return(-tpi/x);
00212 }
00213 z = x*x;
00214 u =
U0[0]+z*(
U0[1]+z*(
U0[2]+z*(
U0[3]+z*
U0[4])));
00215 v = one+z*(
V0[0]+z*(
V0[1]+z*(
V0[2]+z*(
V0[3]+z*
V0[4]))));
00216
return(x*(u/v) + tpi*(
__ieee754_j1(x)*
__ieee754_log(x)-one/x));
00217 }
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
#ifdef __STDC__
00230
static const double pr8[6] = {
00231
#else
00232 static double pr8[6] = {
00233
#endif
00234
0.00000000000000000000e+00,
00235 1.17187499999988647970e-01,
00236 1.32394806593073575129e+01,
00237 4.12051854307378562225e+02,
00238 3.87474538913960532227e+03,
00239 7.91447954031891731574e+03,
00240 };
00241
#ifdef __STDC__
00242
static const double ps8[5] = {
00243
#else
00244 static double ps8[5] = {
00245
#endif
00246
1.14207370375678408436e+02,
00247 3.65093083420853463394e+03,
00248 3.69562060269033463555e+04,
00249 9.76027935934950801311e+04,
00250 3.08042720627888811578e+04,
00251 };
00252
00253
#ifdef __STDC__
00254
static const double pr5[6] = {
00255
#else
00256 static double pr5[6] = {
00257
#endif
00258
1.31990519556243522749e-11,
00259 1.17187493190614097638e-01,
00260 6.80275127868432871736e+00,
00261 1.08308182990189109773e+02,
00262 5.17636139533199752805e+02,
00263 5.28715201363337541807e+02,
00264 };
00265
#ifdef __STDC__
00266
static const double ps5[5] = {
00267
#else
00268 static double ps5[5] = {
00269
#endif
00270
5.92805987221131331921e+01,
00271 9.91401418733614377743e+02,
00272 5.35326695291487976647e+03,
00273 7.84469031749551231769e+03,
00274 1.50404688810361062679e+03,
00275 };
00276
00277
#ifdef __STDC__
00278
static const double pr3[6] = {
00279
#else
00280 static double pr3[6] = {
00281
#endif
00282
3.02503916137373618024e-09,
00283 1.17186865567253592491e-01,
00284 3.93297750033315640650e+00,
00285 3.51194035591636932736e+01,
00286 9.10550110750781271918e+01,
00287 4.85590685197364919645e+01,
00288 };
00289
#ifdef __STDC__
00290
static const double ps3[5] = {
00291
#else
00292 static double ps3[5] = {
00293
#endif
00294
3.47913095001251519989e+01,
00295 3.36762458747825746741e+02,
00296 1.04687139975775130551e+03,
00297 8.90811346398256432622e+02,
00298 1.03787932439639277504e+02,
00299 };
00300
00301
#ifdef __STDC__
00302
static const double pr2[6] = {
00303
#else
00304 static double pr2[6] = {
00305
#endif
00306
1.07710830106873743082e-07,
00307 1.17176219462683348094e-01,
00308 2.36851496667608785174e+00,
00309 1.22426109148261232917e+01,
00310 1.76939711271687727390e+01,
00311 5.07352312588818499250e+00,
00312 };
00313
#ifdef __STDC__
00314
static const double ps2[5] = {
00315
#else
00316 static double ps2[5] = {
00317
#endif
00318
2.14364859363821409488e+01,
00319 1.25290227168402751090e+02,
00320 2.32276469057162813669e+02,
00321 1.17679373287147100768e+02,
00322 8.36463893371618283368e+00,
00323 };
00324
00325
#ifdef __STDC__
00326
static double pone(
double x)
00327 #
else
00328 static double pone(x)
00329 double x;
00330 #endif
00331 {
00332
#ifdef __STDC__
00333
const double *p = 0,*q = 0;
00334
#else
00335
double *p = 0,*q = 0;
00336
#endif
00337
double z,r,s;
00338 int32_t ix;
00339
GET_HIGH_WORD(ix,x);
00340 ix &= 0x7fffffff;
00341
if(ix>=0x40200000) {p =
pr8; q=
ps8;}
00342
else if(ix>=0x40122E8B){p =
pr5; q=
ps5;}
00343
else if(ix>=0x4006DB6D){p =
pr3; q=
ps3;}
00344
else if(ix>=0x40000000){p =
pr2; q=
ps2;}
00345 z = one/(x*x);
00346 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
00347 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
00348
return one+ r/s;
00349 }
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
#ifdef __STDC__
00363
static const double qr8[6] = {
00364
#else
00365 static double qr8[6] = {
00366
#endif
00367
0.00000000000000000000e+00,
00368 -1.02539062499992714161e-01,
00369 -1.62717534544589987888e+01,
00370 -7.59601722513950107896e+02,
00371 -1.18498066702429587167e+04,
00372 -4.84385124285750353010e+04,
00373 };
00374
#ifdef __STDC__
00375
static const double qs8[6] = {
00376
#else
00377 static double qs8[6] = {
00378
#endif
00379
1.61395369700722909556e+02,
00380 7.82538599923348465381e+03,
00381 1.33875336287249578163e+05,
00382 7.19657723683240939863e+05,
00383 6.66601232617776375264e+05,
00384 -2.94490264303834643215e+05,
00385 };
00386
00387
#ifdef __STDC__
00388
static const double qr5[6] = {
00389
#else
00390 static double qr5[6] = {
00391
#endif
00392
-2.08979931141764104297e-11,
00393 -1.02539050241375426231e-01,
00394 -8.05644828123936029840e+00,
00395 -1.83669607474888380239e+02,
00396 -1.37319376065508163265e+03,
00397 -2.61244440453215656817e+03,
00398 };
00399
#ifdef __STDC__
00400
static const double qs5[6] = {
00401
#else
00402 static double qs5[6] = {
00403
#endif
00404
8.12765501384335777857e+01,
00405 1.99179873460485964642e+03,
00406 1.74684851924908907677e+04,
00407 4.98514270910352279316e+04,
00408 2.79480751638918118260e+04,
00409 -4.71918354795128470869e+03,
00410 };
00411
00412
#ifdef __STDC__
00413
static const double qr3[6] = {
00414
#else
00415 static double qr3[6] = {
00416
#endif
00417
-5.07831226461766561369e-09,
00418 -1.02537829820837089745e-01,
00419 -4.61011581139473403113e+00,
00420 -5.78472216562783643212e+01,
00421 -2.28244540737631695038e+02,
00422 -2.19210128478909325622e+02,
00423 };
00424
#ifdef __STDC__
00425
static const double qs3[6] = {
00426
#else
00427 static double qs3[6] = {
00428
#endif
00429
4.76651550323729509273e+01,
00430 6.73865112676699709482e+02,
00431 3.38015286679526343505e+03,
00432 5.54772909720722782367e+03,
00433 1.90311919338810798763e+03,
00434 -1.35201191444307340817e+02,
00435 };
00436
00437
#ifdef __STDC__
00438
static const double qr2[6] = {
00439
#else
00440 static double qr2[6] = {
00441
#endif
00442
-1.78381727510958865572e-07,
00443 -1.02517042607985553460e-01,
00444 -2.75220568278187460720e+00,
00445 -1.96636162643703720221e+01,
00446 -4.23253133372830490089e+01,
00447 -2.13719211703704061733e+01,
00448 };
00449
#ifdef __STDC__
00450
static const double qs2[6] = {
00451
#else
00452 static double qs2[6] = {
00453
#endif
00454
2.95333629060523854548e+01,
00455 2.52981549982190529136e+02,
00456 7.57502834868645436472e+02,
00457 7.39393205320467245656e+02,
00458 1.55949003336666123687e+02,
00459 -4.95949898822628210127e+00,
00460 };
00461
00462
#ifdef __STDC__
00463
static double qone(
double x)
00464 #
else
00465 static double qone(x)
00466 double x;
00467 #endif
00468 {
00469
#ifdef __STDC__
00470
const double *p = 0,*q = 0;
00471
#else
00472
double *p = 0,*q = 0;
00473
#endif
00474
double s,r,z;
00475 int32_t ix;
00476
GET_HIGH_WORD(ix,x);
00477 ix &= 0x7fffffff;
00478
if(ix>=0x40200000) {p =
qr8; q=
qs8;}
00479
else if(ix>=0x40122E8B){p =
qr5; q=
qs5;}
00480
else if(ix>=0x4006DB6D){p =
qr3; q=
qs3;}
00481
else if(ix>=0x40000000){p =
qr2; q=
qs2;}
00482 z = one/(x*x);
00483 r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
00484 s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
00485
return (.375 + r/s)/x;
00486 }