Date: Sun, 12 Feb 95 18:49:29 +0000
This is a transcription of algorithm #5 of Collected Algorithms from
ACM. I took a few liberties to express it in modern ASCII codes, namely,
the 'not equal' construct is expressed as =/= and the exponentiation
is expressed as a ^ character. Hence ((X/2) ^ (n+2*s)) means
n+2*s
X / 2
multiplication operators have been substituted by an * as in most modern
languages.
Jose R. Valverde
European Bioinformatics Institute
txomsy@ebi.ac.uk
------------------------------------------------------------------------
5. BESSEL FUNCTION I, SERIES EXPANSION
Dorothea S. Clarke
General Electric Co., FPLD, Cincinnati 15, Ohio
comment Compute the Bessel function In(X) when n and X
are within the bounds of the series expansion.
The procedure calling statement gives n, X and an
absolute tolerance delta for determining the point at
which the terms of the summation become insig-
nificant. Special case: I0(0) = 1;
procedure I(n, X, delta) =:(Is)
begin
I: s := 0 ; sum := 0
if (n =/= 0) ; go to STRT
if (X = 0) ; begin Is := 1 ; return end
summ := 1 ; go to SURE
STRT: sfac := 1
if (s = 0) ; go to HRE
for t := 1 (1) s
sfac := sfac * t
HRE: snfac := sfac
for t := s + 1 (1) s + n
snfac := snfac * t
summ := sum + ((X/2) ^ (n+2*s)) / sfac * snfac)
SURE: if (delta < abs(summ - sum))
begin s := s + 1 ; sum := summ ; go to STRT end
Is := summ ; return
end