(defun qroot (p n b c eps &key (nmax 20) (itmax 20) (tiny 1.0E-6))
 (declare (type (simple-array double-float (*)) p)) (declare (type fixnum n))
 (declare (type double-float b)) (declare (type double-float c))
 (declare (type double-float eps)) (declare (type fixnum nmax))
 (declare (type fixnum itmax)) (declare (type double-float tiny))
 (prog
  ((q (make-array '(nmax) :element-type 'double-float))
   (d (make-array '(3) :element-type 'double-float))
   (rem (make-array '(nmax) :element-type 'double-float))
   (qq (make-array '(nmax) :element-type 'double-float)) (delc 0.0d0)
   (delb 0.0d0) (div 0.0d0) (rb 0.0d0) (sb 0.0d0) (i 0) (rc 0.0d0) (sc 0.0d0)
   (r 0.0d0) (s 0.0d0) (iter 0)
  )
  (declare (type (simple-array double-float (*)) q))
  (declare (type (simple-array double-float (*)) d))
  (declare (type (simple-array double-float (*)) rem))
  (declare (type (simple-array double-float (*)) qq))
  (declare (type double-float delc)) (declare (type double-float delb))
  (declare (type double-float div)) (declare (type double-float rb))
  (declare (type double-float sb)) (declare (type fixnum i))
  (declare (type double-float rc)) (declare (type double-float sc))
  (declare (type double-float r)) (declare (type double-float s))
  (declare (type fixnum iter)) (fset (fref d 3) 1.0)
  (fdo ((iter 1 (+ iter 1))) ((> iter itmax) nil)
   (tagbody (fset (fref d 2) b) (fset (fref d 1) c)
    (multiple-value-setq (p n d dummy_var q rem) (poldiv p n d 3 q rem))
    (setf s (fref rem 1)) (setf r (fref rem 2))
    (multiple-value-setq (q dummy_var d dummy_var qq rem)
     (poldiv q (+ n (- 1)) d 3 qq rem)
    )
    (setf sc (- (fref rem 1))) (setf rc (- (fref rem 2)))
    (fdo ((i (+ n (- 1)) (+ i (- 1)))) ((> i 1) nil)
     (tagbody (fset (fref q (+ i 1)) (fref q i)))
    )
    (fset (fref q 1) 0.0)
    (multiple-value-setq (q n d dummy_var qq rem) (poldiv q n d 3 qq rem))
    (setf sb (- (fref rem 1))) (setf rb (- (fref rem 2)))
    (setf div (/ 1.0 (+ (* sb rc) (* (* -1 sc) rb))))
    (setf delb (* (+ (* r sc) (* (* -1 s) rc)) div))
    (setf delc (* (+ (* (* -1 r) sb) (* s rb)) div)) (setf b (+ b delb))
    (setf c (+ c delc))
    (if
     (and (or (<= (abs delb) (* eps (abs b))) (< (abs b) tiny))
      (or (<= (abs delc) (* eps (abs c))) (< (abs c) tiny))
     )
     (go end_label)
  )))
  (return (values p n b c eps))
))