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.8 1999/07/02 15:37:44 simonb 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 #ifdef __STDC__ 48 double jn(int n, double x) /* wrapper jn */ 49 #else 50 double jn(n,x) /* wrapper jn */ 51 double x; int n; 52 #endif 53 { 54 #ifdef _IEEE_LIBM 55 return __ieee754_jn(n,x); 56 #else 57 double z; 58 z = __ieee754_jn(n,x); 59 if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; 60 if(fabs(x)>X_TLOSS) { 61 return __kernel_standard((double)n,x,38); /* jn(|x|>X_TLOSS,n) */ 62 } else 63 return z; 64 #endif 65 } 66 67 #ifdef __STDC__ 68 double yn(int n, double x) /* wrapper yn */ 69 #else 70 double yn(n,x) /* wrapper yn */ 71 double x; int n; 72 #endif 73 { 74 #ifdef _IEEE_LIBM 75 return __ieee754_yn(n,x); 76 #else 77 double z; 78 z = __ieee754_yn(n,x); 79 if(_LIB_VERSION == _IEEE_ || isnan(x) ) return z; 80 if(x <= 0.0){ 81 if(x==0.0) 82 /* d= -one/(x-x); */ 83 return __kernel_standard((double)n,x,12); 84 else 85 /* d = zero/(x-x); */ 86 return __kernel_standard((double)n,x,13); 87 } 88 if(x>X_TLOSS) { 89 return __kernel_standard((double)n,x,39); /* yn(x>X_TLOSS,n) */ 90 } else 91 return z; 92 #endif 93 } 94