Function: math-poly-laguerre-root

math-poly-laguerre-root is a byte-compiled function defined in calcalg2.el.gz.

Signature

(math-poly-laguerre-root P X POLISH)

Source Code

;; Defined in /usr/src/emacs/lisp/calc/calcalg2.el.gz
;;; The following routine is from Numerical Recipes, section 9.5.
(defun math-poly-laguerre-root (p x polish)
  (let* ((calc-prefer-frac nil)
	 (calc-symbolic-mode nil)
	 (iters 0)
	 (m (1- (length p)))
	 (try-newt (not polish))
	 ;; (tried-newt nil)
	 b d f x1 dx dxold)
    (while
	(and (or (< (setq iters (1+ iters)) 50)
		 (math-reject-arg x "*Laguerre's method failed to converge"))
	     (let ((err (math-abs-approx (car p)))
		   (abx (math-abs-approx x))
		   (pp p))
	       (setq b (car p)
		     d 0 f 0)
	       (while (setq pp (cdr pp))
		 (setq f (math-add (math-mul x f) d)
		       d (math-add (math-mul x d) b)
		       b (math-add (math-mul x b) (car pp))
		       err (math-add (math-abs-approx b) (math-mul abx err))))
	       (math-lessp (calcFunc-scf err (- -2 calc-internal-prec))
			   (math-abs-approx b)))
	     (or (not (math-zerop d))
		 (not (math-zerop f))
		 (progn
		   (setq x (math-pow (math-neg b) (list 'frac 1 m)))
		   nil))
	     (let* ((g (math-div d b))
		    (g2 (math-sqr g))
		    (h (math-sub g2 (math-mul 2 (math-div f b))))
		    (sq (math-sqrt
			 (math-mul (1- m) (math-sub (math-mul m h) g2))))
		    (gp (math-add g sq))
		    (gm (math-sub g sq)))
	       (if (math-lessp (calcFunc-abssqr gp) (calcFunc-abssqr gm))
		   (setq gp gm))
	       (setq dx (math-div m gp)
		     x1 (math-sub x dx))
	       (if (and try-newt
			(math-lessp (math-abs-approx dx)
				    (calcFunc-scf (math-abs-approx x) -3)))
		   (let ((newt (math-poly-newton-root p x1 7)))
		     (setq ;; tried-newt t
			   try-newt nil)
		     (if (math-zerop (cdr newt))
			 (setq x (car newt) x1 x)
		       (if (math-lessp (cdr newt) '(float 1 -6))
			   (let ((newt2 (math-poly-newton-root
					 p (car newt) 20)))
			     (if (math-zerop (cdr newt2))
				 (setq x (car newt2) x1 x)
			       (setq x (car newt))))))))
	       (not (or (eq x x1)
			(math-nearly-equal x x1))))
	     (let ((cdx (math-abs-approx dx)))
	       (setq x x1
		     ;; tried-newt nil
		     )
	       (prog1
		   (or (<= iters 6)
		       (math-lessp cdx dxold)
		       (progn
			 (if polish
			     (let ((digs (calcFunc-xpon
					  (math-div (math-abs-approx x) cdx))))
			       (calc-record-why
				"*Could not attain full precision")
			       (if (natnump digs)
				   (let ((calc-internal-prec (max 3 digs)))
				     (setq x (math-normalize x))))))
			 nil))
		 (setq dxold cdx)))
	     (or polish
		 (math-lessp (calcFunc-scf (math-abs-approx x)
					   (- calc-internal-prec))
			     dxold))))
    (or (and (math-floatp x)
	     (math-poly-integer-root x))
	x)))