1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
|
/*
floating point Bessel's function of
the first and second kinds and of
integer order.
int n;
double x;
jn(n,x);
returns the value of Jn(x) for all
integer values of n and all real values
of x.
There are no error returns.
Calls j0l, j1l.
For n=0, j0l(x) is called,
for n=1, j1l(x) is called,
for n<x, forward recursion us used starting
from values of j0l(x) and j1l(x).
for n>x, a continued fraction approximation to
j(n,x)/j(n-1,x) is evaluated and then backward
recursion is used starting from a supposed value
for j(n,x). The resulting value of j(0,x) is
compared with the actual value to correct the
supposed value of j(n,x).
yn(n,x) is similar in all respects, except
that forward recursion is used for all
values of n>1.
*/
#include <math.h>
/* long double j0l(), j1l(), y0l(), y1l(); */
long double jnl(int n, long double x)
{
int i, sign;
long double a, b, temp;
long double xsq, t;
/*
n
J (x) = (-1) J (x)
-n n
n
J (-x) = (-1) J (x)
n n
*/
sign = 1;
if(n<0){
n = -n;
if( n & 1 )
sign = -1;
}
if( x < 0.0L )
{
if( n & 1 )
sign = -sign;
x = -x;
}
if(n==0) return(j0l(x));
if(n==1) return(sign*j1l(x));
if(x == 0.L) return(0.L);
if(n>x) goto recurs;
a = j0l(x);
b = j1l(x);
for(i=1;i<n;i++){
temp = b;
b = (2.L*i/x)*b - a;
a = temp;
}
return(sign*b);
recurs:
xsq = x*x;
for(t=0,i=n+16;i>n;i--){
t = xsq/(2.L*i - t);
}
t = x/(2.L*n-t);
a = t;
b = 1;
for(i=n-1;i>0;i--){
temp = b;
b = (2.L*i/x)*b - a;
a = temp;
}
return(sign*t*j0l(x)/b);
}
long double ynl(int n, long double x)
{
int i;
int sign;
long double a, b, temp;
if (x <= 0) {
return(-1.189731495357231765021263853E4932L);
}
sign = 1;
if(n<0){
n = -n;
if(n%2 == 1) sign = -1;
}
if(n==0) return(y0l(x));
if(n==1) return(sign*y1l(x));
a = y0l(x);
b = y1l(x);
for(i=1;i<n;i++){
temp = b;
b = (2.L*i/x)*b - a;
a = temp;
}
return(sign*b);
}
|
