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 #include <sys/cdefs.h> 14 #if defined(LIBM_SCCS) && !defined(lint) 15 __RCSID("$NetBSD: w_jn.c,v 1.9 2002/05/26 22:02:01 wiz Exp $"); 16 #endif 17 18 /* 19 * wrapper jn(int n, double x), yn(int n, double x) 20 * floating point Bessel's function of the 1st and 2nd kind 21 * of order n 22 * 23 * Special cases: 24 * y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal; 25 * y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal. 26 * Note 2. About jn(n,x), yn(n,x) 27 * For n=0, j0(x) is called, 28 * for n=1, j1(x) is called, 29 * for n<x, forward recursion us used starting 30 * from values of j0(x) and j1(x). 31 * for n>x, a continued fraction approximation to 32 * j(n,x)/j(n-1,x) is evaluated and then backward 33 * recursion is used starting from a supposed value 34 * for j(n,x). The resulting value of j(0,x) is 35 * compared with the actual value to correct the 36 * supposed value of j(n,x). 37 * 38 * yn(n,x) is similar in all respects, except 39 * that forward recursion is used for all 40 * values of n>1. 41 * 42 */ 43 44 #include "math.h" 45 #include "math_private.h" 46 47 double 48 jn(int n, double x) /* wrapper jn */ 49 { 50 #ifdef _IEEE_LIBM 51 return __ieee754_jn(n,x); 52 #else 53 double z; 54 z = __ieee754_jn(n,x); 55 if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; 56 if(fabs(x)>X_TLOSS) { 57 return __kernel_standard((double)n,x,38); /* jn(|x|>X_TLOSS,n) */ 58 } else 59 return z; 60 #endif 61 } 62 63 double 64 yn(int n, double x) /* wrapper yn */ 65 { 66 #ifdef _IEEE_LIBM 67 return __ieee754_yn(n,x); 68 #else 69 double z; 70 z = __ieee754_yn(n,x); 71 if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; 72 if(x <= 0.0){ 73 if(x==0.0) 74 /* d= -one/(x-x); */ 75 return __kernel_standard((double)n,x,12); 76 else 77 /* d = zero/(x-x); */ 78 return __kernel_standard((double)n,x,13); 79 } 80 if(x>X_TLOSS) { 81 return __kernel_standard((double)n,x,39); /* yn(x>X_TLOSS,n) */ 82 } else 83 return z; 84 #endif 85 } 86