1 2 /* @(#)w_jn.c 1.3 95/01/18 */ 3 /* 4 * ==================================================== 5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6 * 7 * Developed at SunSoft, a Sun Microsystems, Inc. business. 8 * Permission to use, copy, modify, and distribute this 9 * software is freely granted, provided that this notice 10 * is preserved. 11 * ==================================================== 12 */ 13 14 /* 15 * wrapper ieee_jn(int n, double x), ieee_yn(int n, double x) 16 * floating point Bessel's function of the 1st and 2nd kind 17 * of order n 18 * 19 * Special cases: 20 * y0(0)=ieee_y1(0)=ieee_yn(n,0) = -inf with division by zero signal; 21 * y0(-ve)=ieee_y1(-ve)=ieee_yn(n,-ve) are NaN with invalid signal. 22 * Note 2. About ieee_jn(n,x), ieee_yn(n,x) 23 * For n=0, ieee_j0(x) is called, 24 * for n=1, ieee_j1(x) is called, 25 * for n<x, forward recursion us used starting 26 * from values of ieee_j0(x) and ieee_j1(x). 27 * for n>x, a continued fraction approximation to 28 * j(n,x)/j(n-1,x) is evaluated and then backward 29 * recursion is used starting from a supposed value 30 * for j(n,x). The resulting value of j(0,x) is 31 * compared with the actual value to correct the 32 * supposed value of j(n,x). 33 * 34 * yn(n,x) is similar in all respects, except 35 * that forward recursion is used for all 36 * values of n>1. 37 * 38 */ 39 40 #include "fdlibm.h" 41 42 #ifdef __STDC__ ieee_jn(int n,double x)43 double ieee_jn(int n, double x) /* wrapper jn */ 44 #else 45 double ieee_jn(n,x) /* wrapper jn */ 46 double x; int n; 47 #endif 48 { 49 #ifdef _IEEE_LIBM 50 return __ieee754_jn(n,x); 51 #else 52 double z; 53 z = __ieee754_jn(n,x); 54 if(_LIB_VERSION == _IEEE_ || ieee_isnan(x) ) return z; 55 if(ieee_fabs(x)>X_TLOSS) { 56 return __kernel_standard((double)n,x,38); /* ieee_jn(|x|>X_TLOSS,n) */ 57 } else 58 return z; 59 #endif 60 } 61 62 #ifdef __STDC__ ieee_yn(int n,double x)63 double ieee_yn(int n, double x) /* wrapper yn */ 64 #else 65 double ieee_yn(n,x) /* wrapper yn */ 66 double x; int n; 67 #endif 68 { 69 #ifdef _IEEE_LIBM 70 return __ieee754_yn(n,x); 71 #else 72 double z; 73 z = __ieee754_yn(n,x); 74 if(_LIB_VERSION == _IEEE_ || ieee_isnan(x) ) return z; 75 if(x <= 0.0){ 76 if(x==0.0) 77 /* d= -one/(x-x); */ 78 return __kernel_standard((double)n,x,12); 79 else 80 /* d = zero/(x-x); */ 81 return __kernel_standard((double)n,x,13); 82 } 83 if(x>X_TLOSS) { 84 return __kernel_standard((double)n,x,39); /* ieee_yn(x>X_TLOSS,n) */ 85 } else 86 return z; 87 #endif 88 } 89