;;;-*- mode: lisp; package: ccpl -*- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; File: nr.lisp
;;
;;
;; Date of Creation: 1998.08.18
;;
;;
;; Purpose: Parts of Numerical Recipies
;;
;; Version: 0.00
;;
;; History of Changes:
;;;
;;; translated by Bill Schottstaedt <bil@ccrma.Stanford.EDU> from "Numerical Recipes
;;; in Pascal" by Press, Flannery, Teukolsky, Vetterling
;;; Cambridge Univ Press, 1989. All functions and variables are named the same.
;;;; ------------------------------------------------------------------------------------------
;;;; Root Finding
(eval-when (:compile-toplevel :load-toplevel :execute)
(in-package :ccpl))
(declaim (optimize (safty 0)(speed 3)(space 3)(debug 0)))
(defun sqr (x) (* x x))
(defun laguer (a m x eps polish)
(declare (inline sqr))
(let ((epss 6.0e-8)
(mixt 100)
(dxold (abs x)))
(loop for iter from 1 to mixt do
(let* ((b (aref a m))
(err (abs b))
(d #C(0 0))
(f #C(0 0))
(g #C(0 0))
(g2 #C(0 0))
(gp #C(0 0))
(gm #C(0 0))
(h #C(0 0))
(cdum #C(0 0))
(sq #C(0 0))
(dx #C(0 0))
(cdx #C(0 0))
(x1 #C(0 0))
(abx (abs x)))
(loop for j from (1- m) downto 0 do
(let ((dum (realpart f)))
(setf f (complex (+ (* (realpart x) (realpart f))
(realpart d)
(- (* (imagpart x) (imagpart f))))
(+ (* (realpart x) (imagpart f))
(imagpart d)
(* (imagpart x) dum))))
(setf dum (realpart d))
(setf d (complex (+ (* (realpart x) (realpart d))
(realpart b)
(- (* (imagpart x) (imagpart d))))
(+ (* (realpart x) (imagpart d))
(imagpart b)
(* (imagpart x) dum))))
(setf dum (realpart b))
(setf b (complex (+ (* (realpart x) (realpart b))
(realpart (aref a j))
(- (* (imagpart x) (imagpart b))))
(+ (* (realpart x) (imagpart b))
(imagpart (aref a j))
(* (imagpart x) dum))))
(setf err (+ (abs b) (* abx err)))))
(setf err (* epss err))
(if (<= (abs b) err) (return-from laguer x))
(setf g (/ d b))
(setf g2 (complex (- (sqr (realpart g)) (sqr (imagpart g)))
(* 2.0 (realpart g) (imagpart g))))
(setf cdum (/ f b))
(setf h (complex (- (realpart g2) (* 2.0 (realpart cdum)))
(- (imagpart g2) (* 2.0 (imagpart cdum)))))
(setf cdum (complex (* (1- m) (- (* m (realpart h)) (realpart g2)))
(* (1- m) (- (* m (imagpart h)) (imagpart g2)))))
(setf sq (sqrt cdum))
(setf gp (+ g sq))
(setf gm (- g sq))
(if (< (abs gp) (abs gm)) (setf gp gm))
(setf cdum (complex m 0))
(setf dx (/ cdum gp))
(setf x1 (- x dx))
(if (= x x1) (return-from laguer x))
(setf x x1)
(setf cdx (abs dx))
(setf dxold cdx)
(if (and (not polish)
(<= cdx (* eps (abs x))))
(return-from laguer x))))))
(defun zroots (a m roots polish)
"Find root of polymonic equations
Given the degree m and the m+1 complex coefficients a[0..m] of the polynoial
this routine successivley calles laguer and finds all m complex roots i roots[1..m] The
logical variable polich should be input as T if polishing is desired, NIL if the roots
will be polished by other means."
(let ((eps 2.0e-12)
(ad (make-array (+ m 2) :element-type 'complex))
(b #C(0 0))
(cc #C(0 0))
(dum 0.0))
(declare (type complex a cc)
(type (vector complex) ad)
(dynamic-extent ad))
(loop for j from 0 to m do
(setf (aref ad j) (aref a j)))
(loop for j from m downto 1 do
(let ((x (laguer ad j #C(0 0) eps nil)))
(if (<= (abs (imagpart x))
(* 2.0 eps (abs (realpart x))))
(setf x (complex (realpart x) 0.0)))
(setf (aref roots j) x)
(setf b (aref ad j))
(loop for jj from (1- j) downto 0 do
(setf cc (aref ad jj))
(setf (aref ad jj) b)
(setf dum (realpart b))
(setf b (complex (+ (* (realpart b) (realpart x))
(realpart cc)
(- (* (imagpart b) (imagpart x))))
(+ (* dum (imagpart x))
(imagpart cc)
(* (imagpart b) (realpart x))))))))
(if polish
(loop for j from 1 to m do
(setf (aref roots j) (laguer a m (aref roots j) eps t))))
(loop for j from 2 to m do
(let ((x (aref roots j))
(i (1- j)))
(loop while (and (>= i 1)
(> (realpart (aref roots i)) (realpart x))) do
(setf (aref roots (1+ i)) (aref roots i))
(decf i))
(setf (aref roots (1+ i)) x)))))
