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_jn.c,v 1.9 1995/05/10 20:45:34 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
#include "math.h"
00044
#include "mathP.h"
00045
00046
#ifdef __STDC__
00047
static const double
00048
#else
00049
static double
00050
#endif
00051 invsqrtpi= 5.64189583547756279280e-01,
00052 two = 2.00000000000000000000e+00,
00053 one = 1.00000000000000000000e+00;
00054
00055
#ifdef __STDC__
00056
static const double zero = 0.00000000000000000000e+00;
00057
#else
00058 static double zero = 0.00000000000000000000e+00;
00059
#endif
00060
00061
#ifdef __STDC__
00062
double __ieee754_jn(
int n,
double x)
00063 #
else
00064 double __ieee754_jn(n,x)
00065 int n;
double x;
00066 #endif
00067 {
00068 int32_t i,hx,ix,lx, sgn;
00069
double a, b, temp, di;
00070
double z, w;
00071
00072
00073
00074
00075
EXTRACT_WORDS(hx,lx,x);
00076 ix = 0x7fffffff&hx;
00077
00078
if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000)
return x+x;
00079
if(n<0){
00080 n = -n;
00081 x = -x;
00082 hx ^= 0x80000000;
00083 }
00084
if(n==0)
return(
__ieee754_j0(x));
00085
if(n==1)
return(
__ieee754_j1(x));
00086 sgn = (n&1)&(hx>>31);
00087 x =
fabs(x);
00088
if((ix|lx)==0||ix>=0x7ff00000)
00089 b =
zero;
00090
else if((
double)n<=x) {
00091
00092
if(ix>=0x52D00000) {
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
switch(n&3) {
00107
case 0: temp =
cos(x)+
sin(x);
break;
00108
case 1: temp = -
cos(x)+
sin(x);
break;
00109
case 2: temp = -
cos(x)-
sin(x);
break;
00110
case 3: temp =
cos(x)-
sin(x);
break;
00111
default: temp = 0.0;
00112 }
00113 b =
invsqrtpi*temp/
sqrt(x);
00114 }
else {
00115 a =
__ieee754_j0(x);
00116 b =
__ieee754_j1(x);
00117
for(i=1;i<n;i++){
00118 temp = b;
00119 b = b*((
double)(i+i)/x) - a;
00120 a = temp;
00121 }
00122 }
00123 }
else {
00124
if(ix<0x3e100000) {
00125
00126
00127
00128
if(n>33)
00129 b =
zero;
00130
else {
00131 temp = x*0.5; b = temp;
00132
for (a=one,i=2;i<=n;i++) {
00133 a *= (
double)i;
00134 b *= temp;
00135 }
00136 b = b/a;
00137 }
00138 }
else {
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
double t,v;
00169
double q0,q1,h,tmp; int32_t k,m;
00170 w = (n+n)/(
double)x; h = 2.0/(
double)x;
00171 q0 = w; z = w+h; q1 = w*z - 1.0; k=1;
00172
while(q1<1.0e9) {
00173 k += 1; z += h;
00174 tmp = z*q1 - q0;
00175 q0 = q1;
00176 q1 = tmp;
00177 }
00178 m = n+n;
00179
for(t=
zero, i = 2*(n+k); i>=m; i -= 2) t = one/(i/x-t);
00180 a = t;
00181 b = one;
00182
00183
00184
00185
00186
00187
00188
00189
00190 tmp = n;
00191 v = two/x;
00192 tmp = tmp*
__ieee754_log(
fabs(v*tmp));
00193
if(tmp<7.09782712893383973096e+02) {
00194
for(i=n-1,di=(
double)(i+i);i>0;i--){
00195 temp = b;
00196 b *= di;
00197 b = b/x - a;
00198 a = temp;
00199 di -= two;
00200 }
00201 }
else {
00202
for(i=n-1,di=(
double)(i+i);i>0;i--){
00203 temp = b;
00204 b *= di;
00205 b = b/x - a;
00206 a = temp;
00207 di -= two;
00208
00209
if(b>1e100) {
00210 a /= b;
00211 t /= b;
00212 b = one;
00213 }
00214 }
00215 }
00216 b = (t*
__ieee754_j0(x)/b);
00217 }
00218 }
00219
if(sgn==1)
return -b;
else return b;
00220 }
00221
00222
#ifdef __STDC__
00223
double __ieee754_yn(
int n,
double x)
00224 #
else
00225 double __ieee754_yn(n,x)
00226 int n;
double x;
00227 #endif
00228 {
00229 int32_t i,hx,ix,lx;
00230 int32_t sign;
00231
double a, b, temp;
00232
00233
EXTRACT_WORDS(hx,lx,x);
00234 ix = 0x7fffffff&hx;
00235
00236
if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000)
return x+x;
00237
if((ix|lx)==0)
return -one/
zero;
00238
if(hx<0)
return zero/
zero;
00239 sign = 1;
00240
if(n<0){
00241 n = -n;
00242 sign = 1 - ((n&1)<<1);
00243 }
00244
if(n==0)
return(
__ieee754_y0(x));
00245
if(n==1)
return(sign*
__ieee754_y1(x));
00246
if(ix==0x7ff00000)
return zero;
00247
if(ix>=0x52D00000) {
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
switch(n&3) {
00262
case 0: temp =
sin(x)-
cos(x);
break;
00263
case 1: temp = -
sin(x)-
cos(x);
break;
00264
case 2: temp = -
sin(x)+
cos(x);
break;
00265
case 3: temp =
sin(x)+
cos(x);
break;
00266
default: temp = 0.0;
00267 }
00268 b =
invsqrtpi*temp/
sqrt(x);
00269 }
else {
00270 u_int32_t high;
00271 a =
__ieee754_y0(x);
00272 b =
__ieee754_y1(x);
00273
00274
GET_HIGH_WORD(high,b);
00275
for(i=1;i<n&&high!=0xfff00000;i++){
00276 temp = b;
00277 b = ((
double)(i+i)/x)*b - a;
00278
GET_HIGH_WORD(high,b);
00279 a = temp;
00280 }
00281 }
00282
if(sign>0)
return b;
else return -b;
00283 }