1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
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))))))))
|
