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