+;;; bug #389 (Rick Taube sbcl-devel)
+(defun bes-jn (unn ux)
+ (let ((nn unn) (x ux))
+ (let* ((n (floor (abs nn)))
+ (besn
+ (if (= n 0)
+ (bes-j0 x)
+ (if (= n 1)
+ (bes-j1 x)
+ (if (zerop x)
+ 0.0
+ (let ((iacc 40)
+ (ans 0.0)
+ (bigno 1.0e+10)
+ (bigni 1.0e-10))
+ (if (> (abs x) n)
+ (do ((tox (/ 2.0 (abs x)))
+ (bjm (bes-j0 (abs x)))
+ (bj (bes-j1 (abs x)))
+ (j 1 (+ j 1))
+ (bjp 0.0))
+ ((= j n) (setf ans bj))
+ (setf bjp (- (* j tox bj) bjm))
+ (setf bjm bj)
+ (setf bj bjp))
+ (let ((tox (/ 2.0 (abs x)))
+ (m
+ (* 2
+ (floor
+ (/ (+ n (sqrt (* iacc n)))
+ 2))))
+ (jsum 0.0)
+ (bjm 0.0)
+ (sum 0.0)
+ (bjp 0.0)
+ (bj 1.0))
+ (do ((j m (- j 1)))
+ ((= j 0))
+ (setf bjm (- (* j tox bj) bjp))
+ (setf bjp bj)
+ (setf bj bjm)
+ (when (> (abs bj) bigno)
+ (setf bj (* bj bigni))
+ (setf bjp (* bjp bigni))
+ (setf ans (* ans bigni))
+ (setf sum (* sum bigni)))
+ (if (not (= 0 jsum)) (incf sum bj))
+ (setf jsum (- 1 jsum))
+ (if (= j n) (setf ans bjp)))
+ (setf sum (- (* 2.0 sum) bj))
+ (setf ans (/ ans sum))))
+ (if (and (minusp x) (oddp n))
+ (- ans)
+ ans)))))))
+ (if (and (minusp nn) (oddp nn)) (- besn) besn))))
+