summaryrefslogtreecommitdiffstats
path: root/root.scm
diff options
context:
space:
mode:
Diffstat (limited to 'root.scm')
-rw-r--r--root.scm149
1 files changed, 149 insertions, 0 deletions
diff --git a/root.scm b/root.scm
new file mode 100644
index 0000000..5ba78c1
--- /dev/null
+++ b/root.scm
@@ -0,0 +1,149 @@
+;;;"root.scm" Newton's and Laguerre's methods for finding roots.
+;Copyright (C) 1996 Aubrey Jaffer
+;
+;Permission to copy this software, to redistribute it, and to use it
+;for any purpose is granted, subject to the following restrictions and
+;understandings.
+;
+;1. Any copy made of this software must include this copyright notice
+;in full.
+;
+;2. I have made no warrantee or representation that the operation of
+;this software will be error-free, and I am under no obligation to
+;provide any services, by way of maintenance, update, or otherwise.
+;
+;3. In conjunction with products arising from the use of this
+;material, there shall be no use of my name in any advertising,
+;promotional, or sales literature without prior written consent in
+;each case.
+
+;;;; Newton's Method explained in:
+;;; D. E. Knuth, "The Art of Computer Programming", Vol 2 /
+;;; Seminumerical Algorithms, Reading Massachusetts, Addison-Wesley
+;;; Publishing Company, 2nd Edition, p. 510
+
+(define (newton:find-integer-root f df/dx x_0)
+ (let loop ((x x_0) (fx (f x_0)))
+ (cond
+ ((zero? fx) x)
+ (else
+ (let ((df (df/dx x)))
+ (cond
+ ((zero? df) #f) ; stuck at local min/max
+ (else
+ (let* ((delta (quotient (+ fx (quotient df 2)) df))
+ (next-x (cond ((not (zero? delta)) (- x delta))
+ ((positive? fx) (- x 1))
+ (else (- x -1))))
+ (next-fx (f next-x)))
+ (cond ((>= (abs next-fx) (abs fx)) x)
+ (else (loop next-x next-fx)))))))))))
+
+(define (integer-sqrt y)
+ (newton:find-integer-root (lambda (x) (- (* x x) y))
+ (lambda (x) (* 2 x))
+ (ash 1 (quotient (integer-length y) 2))))
+
+(define (newton:find-root f df/dx x_0 prec)
+ (if (and (negative? prec) (integer? prec))
+ (let loop ((x x_0) (fx (f x_0)) (count prec))
+ (cond ((zero? count) x)
+ (else (let ((df (df/dx x)))
+ (cond ((zero? df) #f) ; stuck at local min/max
+ (else (let* ((next-x (- x (/ fx df)))
+ (next-fx (f next-x)))
+ (cond ((= next-x x) x)
+ ((> (abs next-fx) (abs fx)) #f)
+ (else (loop next-x next-fx
+ (+ 1 count)))))))))))
+ (let loop ((x x_0) (fx (f x_0)))
+ (cond ((< (abs fx) prec) x)
+ (else (let ((df (df/dx x)))
+ (cond ((zero? df) #f) ; stuck at local min/max
+ (else (let* ((next-x (- x (/ fx df)))
+ (next-fx (f next-x)))
+ (cond ((= next-x x) x)
+ ((> (abs next-fx) (abs fx)) #f)
+ (else (loop next-x next-fx))))))))))))
+
+;;; H. J. Orchard, "The Laguerre Method for Finding the Zeros of
+;;; Polynomials", IEEE Transactions on Circuits and Systems, Vol. 36,
+;;; No. 11, November 1989, pp 1377-1381.
+
+(define (laguerre:find-root f df/dz ddf/dz^2 z_0 prec)
+ (if (and (negative? prec) (integer? prec))
+ (let loop ((z z_0) (fz (f z_0)) (count prec))
+ (cond ((zero? count) z)
+ (else
+ (let* ((df (df/dz z))
+ (ddf (ddf/dz^2 z))
+ (disc (sqrt (- (* df df) (* fz ddf)))))
+ (if (zero? disc)
+ #f
+ (let* ((next-z
+ (- z (/ fz (if (negative? (+ (* (real-part df)
+ (real-part disc))
+ (* (imag-part df)
+ (imag-part disc))))
+ (- disc) disc))))
+ (next-fz (f next-z)))
+ (cond ((>= (magnitude next-fz) (magnitude fz)) z)
+ (else (loop next-z next-fz (+ 1 count))))))))))
+ (let loop ((z z_0) (fz (f z_0)) (delta-z #f))
+ (cond ((< (magnitude fz) prec) z)
+ (else
+ (let* ((df (df/dz z))
+ (ddf (ddf/dz^2 z))
+ (disc (sqrt (- (* df df) (* fz ddf)))))
+ (print 'disc disc)
+ (if (zero? disc)
+ #f
+ (let* ((next-z
+ (- z (/ fz (if (negative? (+ (* (real-part df)
+ (real-part disc))
+ (* (imag-part df)
+ (imag-part disc))))
+ (- disc) disc))))
+ (next-delta-z (magnitude (- next-z z))))
+ (print 'next-z next-z )
+ (print '(f next-z) (f next-z))
+ (print 'delta-z delta-z 'next-delta-z next-delta-z)
+ (cond ((zero? next-delta-z) z)
+ ((and delta-z (>= next-delta-z delta-z)) z)
+ (else
+ (loop next-z (f next-z) next-delta-z)))))))))))
+
+(define (laguerre:find-polynomial-root deg f df/dz ddf/dz^2 z_0 prec)
+ (if (and (negative? prec) (integer? prec))
+ (let loop ((z z_0) (fz (f z_0)) (count prec))
+ (cond ((zero? count) z)
+ (else
+ (let* ((df (df/dz z))
+ (ddf (ddf/dz^2 z))
+ (tmp (* (+ deg -1) df))
+ (sqrt-H (sqrt (- (* tmp tmp) (* deg (+ deg -1) fz ddf))))
+ (df+sqrt-H (+ df sqrt-H))
+ (df-sqrt-H (- df sqrt-H))
+ (next-z
+ (- z (/ (* deg fz)
+ (if (>= (magnitude df+sqrt-H)
+ (magnitude df-sqrt-H))
+ df+sqrt-H
+ df-sqrt-H)))))
+ (loop next-z (f next-z) (+ 1 count))))))
+ (let loop ((z z_0) (fz (f z_0)))
+ (cond ((< (magnitude fz) prec) z)
+ (else
+ (let* ((df (df/dz z))
+ (ddf (ddf/dz^2 z))
+ (tmp (* (+ deg -1) df))
+ (sqrt-H (sqrt (- (* tmp tmp) (* deg (+ deg -1) fz ddf))))
+ (df+sqrt-H (+ df sqrt-H))
+ (df-sqrt-H (- df sqrt-H))
+ (next-z
+ (- z (/ (* deg fz)
+ (if (>= (magnitude df+sqrt-H)
+ (magnitude df-sqrt-H))
+ df+sqrt-H
+ df-sqrt-H)))))
+ (loop next-z (f next-z))))))))