1 /* md_yn.c
2 *
3 * Bessel function of second kind of integer order
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double x, y, md_yn();
10 * int n;
11 *
12 * y = md_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 * md_y0() and md_y1().
24 *
25 * If n = 0 or 1 the routine for md_y0 or md_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 * md_yn singularity x = 0 MAXNUM
44 * md_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.8: June, 2000
52 Copyright 1984, 1987, 2000 by Stephen L. Moshier
53 */
54
55 #include "mconf.h"
56 #ifdef ANSIPROT
57 extern double md_y0 ( double );
58 extern double md_y1 ( double );
59 extern double md_log ( double );
60 #else
61 double md_y0(), md_y1(), md_log();
62 #endif
63 extern double MAXNUM, MAXLOG;
64
md_yn(n,x)65 double md_yn( n, x )
66 int n;
67 double x;
68 {
69 double an, anm1, anm2, r;
70 int k, sign;
71
72 if( n < 0 )
73 {
74 n = -n;
75 if( (n & 1) == 0 ) /* -1**n */
76 sign = 1;
77 else
78 sign = -1;
79 }
80 else
81 sign = 1;
82
83
84 if( n == 0 )
85 return( sign * md_y0(x) );
86 if( n == 1 )
87 return( sign * md_y1(x) );
88
89 /* test for overflow */
90 if( x <= 0.0 )
91 {
92 mtherr( "md_yn", SING );
93 return( -MAXNUM );
94 }
95
96 /* forward recurrence on n */
97
98 anm2 = md_y0(x);
99 anm1 = md_y1(x);
100 k = 1;
101 r = 2 * k;
102 do
103 {
104 an = r * anm1 / x - anm2;
105 anm2 = anm1;
106 anm1 = an;
107 r += 2.0;
108 ++k;
109 }
110 while( k < n );
111
112
113 return( sign * an );
114 }
115