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 #ifndef lint 14 static char rcsid[] = "$Id: w_jn.c,v 1.3 1994/02/18 02:27:45 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 45 #ifdef __STDC__ 46 double jn(int n, double x) /* wrapper jn */ 47 #else 48 double jn(n,x) /* wrapper jn */ 49 double x; int n; 50 #endif 51 { 52 #ifdef _IEEE_LIBM 53 return __ieee754_jn(n,x); 54 #else 55 double z; 56 z = __ieee754_jn(n,x); 57 if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; 58 if(fabs(x)>X_TLOSS) { 59 return __kernel_standard((double)n,x,38); /* jn(|x|>X_TLOSS,n) */ 60 } else 61 return z; 62 #endif 63 } 64 65 #ifdef __STDC__ 66 double yn(int n, double x) /* wrapper yn */ 67 #else 68 double yn(n,x) /* wrapper yn */ 69 double x; int n; 70 #endif 71 { 72 #ifdef _IEEE_LIBM 73 return __ieee754_yn(n,x); 74 #else 75 double z; 76 z = __ieee754_yn(n,x); 77 if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; 78 if(x <= 0.0){ 79 if(x==0.0) 80 /* d= -one/(x-x); */ 81 return __kernel_standard((double)n,x,12); 82 else 83 /* d = zero/(x-x); */ 84 return __kernel_standard((double)n,x,13); 85 } 86 if(x>X_TLOSS) { 87 return __kernel_standard((double)n,x,39); /* yn(x>X_TLOSS,n) */ 88 } else 89 return z; 90 #endif 91 } 92