(defun bessi (n x &key (iacc 40) (bigno 1.0E10) (bigni 1.0E-10))
 (declare (type fixnum n)) (declare (type double-float x))
 (declare (type fixnum iacc)) (declare (type double-float bigno))
 (declare (type double-float bigni))
 (prog
  ((bessi 0.0d0) (bim 0.0d0) (j 0) (m 0) (bi 0.0d0) (bip 0.0d0) (tox 0.0d0))
  (declare (type double-float bessi)) (declare (type double-float bim))
  (declare (type fixnum j)) (declare (type fixnum m))
  (declare (type double-float bi)) (declare (type double-float bip))
  (declare (type double-float tox))
  (if (< n 2) (error "bad argument N in BESSI")) (setf tox (/ 2.0 x))
  (setf bip 0.0) (setf bi 1.0)
  (setf m (* 2 (+ n (int (sqrt (float (* iacc n)))))))
  (fdo ((j m (+ j (- 1)))) ((> j 1) nil)
   (tagbody (setf bim (+ bip (* (* (float j) tox) bi))) (setf bip bi)
    (setf bi bim)
    (cond
     ((> (abs bi) bigno) (setf bessi (* bessi bigni)) (setf bi (* bi bigni))
      (setf bip (* bip bigni))
    ))
    (if (= j n) (setf bessi bip))
  ))
  (setf bessi (/ (* bessi (bessi0 x)) bi)) (return bessi)
))