1 /* yn.c
2 *
3 * Bessel function of second kind of integer order
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, y, yn();
10 * int n;
11 *
12 * y = yn( n, x );
13 *
14 *
15 *
16 * DESCRIPTION:
17 *
18 * Returns Bessel function of order n, where n is a
19 * (possibly negative) integer.
20 *
21 * The function is evaluated by forward recurrence on
22 * n, starting with values computed by the routines
23 * y0() and y1().
24 *
25 * If n = 0 or 1 the routine for y0 or y1 is called
26 * directly.
27 *
28 *
29 *
30 * ACCURACY:
31 *
32 *
33 * Absolute error, except relative
34 * when y > 1:
35 * arithmetic domain # trials peak rms
36 * DEC 0, 30 2200 2.9e-16 5.3e-17
37 * IEEE 0, 30 30000 3.4e-15 4.3e-16
38 *
39 *
40 * ERROR MESSAGES:
41 *
42 * message condition value returned
43 * yn singularity x = 0 MAXNUM
44 * yn overflow MAXNUM
45 *
46 * Spot checked against tables for x, n between 0 and 100.
47 *
48 */
49
50 /*
51 Cephes Math Library Release 2.1: December, 1988
52 Copyright 1984, 1987 by Stephen L. Moshier
53 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
54 */
55
56 #include "mconf.h"
57 #include "cephes.h"
58
59 #ifndef HAVE_YN
60
61 extern double MAXNUM, MAXLOG;
62
yn(n,x)63 double yn( n, x )
64 int n;
65 double x;
66 {
67 double an, anm1, anm2, r;
68 int k, sign;
69
70 if( n < 0 )
71 {
72 n = -n;
73 if( (n & 1) == 0 ) /* -1**n */
74 sign = 1;
75 else
76 sign = -1;
77 }
78 else
79 sign = 1;
80
81
82 if( n == 0 )
83 return( sign * y0(x) );
84 if( n == 1 )
85 return( sign * y1(x) );
86
87 /* test for overflow */
88 if( x <= 0.0 )
89 {
90 mtherr( "yn", SING );
91 return( -MAXNUM );
92 }
93
94 /* forward recurrence on n */
95
96 anm2 = y0(x);
97 anm1 = y1(x);
98 k = 1;
99 r = 2 * k;
100 do
101 {
102 an = r * anm1 / x - anm2;
103 anm2 = anm1;
104 anm1 = an;
105 r += 2.0;
106 ++k;
107 }
108 while( k < n );
109
110
111 return( sign * an );
112 }
113
114 #endif /* HAVE_YN */
115
116