(defun bessk1 (x) (declare (type double-float x)) (prog ((y 0.0d0) (p1 0.0d0) (p2 0.0d0) (p3 0.0d0) (p4 0.0d0) (p5 0.0d0) (p6 0.0d0) (p7 0.0d0) (q1 0.0d0) (q2 0.0d0) (q3 0.0d0) (q4 0.0d0) (q5 0.0d0) (q6 0.0d0) (q7 0.0d0) (bessk1 0.0d0) ) (declare (type float y)) (declare (type float p1)) (declare (type float p2)) (declare (type float p3)) (declare (type float p4)) (declare (type float p5)) (declare (type float p6)) (declare (type float p7)) (declare (type float q1)) (declare (type float q2)) (declare (type float q3)) (declare (type float q4)) (declare (type float q5)) (declare (type float q6)) (declare (type float q7)) (declare (type double-float bessk1)) (setq p1 1.0d0) (setq q1 1.25331414d0) (cond ((<= x 2.0) (setf y (/ (* x x) 4.0)) (setf bessk1 (+ (* (log (/ x 2.0)) (bessi1 x)) (* (/ 1.0 x) (+ p1 (* y (+ p2 (* y (+ p3 (* y (+ p4 (* y (+ p5 (* y (+ p6 (* y p7)))))))))) )))))) (t (setf y (/ 2.0 x)) (setf bessk1 (* (/ (exp (- x)) (sqrt x)) (+ q1 (* y (+ q2 (* y (+ q3 (* y (+ q4 (* y (+ q5 (* y (+ q6 (* y q7)))))))))) )))))) (return bessk1) ))