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