# Subversion RepositoriesKolibri OS

Rev
1. /* Copyright (C) 1994 DJ Delorie, see COPYING.DJ for details */
2. /* @(#)w_jn.c 5.1 93/09/24 */
3. /*
4.  * ====================================================
6.  *
7.  * Developed at SunPro, a Sun Microsystems, Inc. business.
8.  * Permission to use, copy, modify, and distribute this
9.  * software is freely granted, provided that this notice
10.  * is preserved.
11.  * ====================================================
12.  */
13.
14. #if defined(LIBM_SCCS) && !defined(lint)
15. static char rcsid[] = "\$Id: w_jn.c,v 1.4 1994/08/10 20:34:42 jtc 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.