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