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 47double 48jn(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 63double 64yn(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