module/slib/root.scm

   1 ;;;"root.scm" Newton's and Laguerre's methods for finding roots.
   2 ;Copyright (C) 1996, 1997 Aubrey Jaffer
   3 ;
   4 ;Permission to copy this software, to redistribute it, and to use it
   5 ;for any purpose is granted, subject to the following restrictions and
   6 ;understandings.
   7 ;
   8 ;1.  Any copy made of this software must include this copyright notice
   9 ;in full.
  10 ;
  11 ;2.  I have made no warrantee or representation that the operation of
  12 ;this software will be error-free, and I am under no obligation to
  13 ;provide any services, by way of maintenance, update, or otherwise.
  14 ;
  15 ;3.  In conjunction with products arising from the use of this
  16 ;material, there shall be no use of my name in any advertising,
  17 ;promotional, or sales literature without prior written consent in
  18 ;each case.
  19
  20 (require 'logical)
  21
  22 ;;;; Newton's Method explained in:
  23 ;;; D. E. Knuth, "The Art of Computer Programming", Vol 2 /
  24 ;;; Seminumerical Algorithms, Reading Massachusetts, Addison-Wesley
  25 ;;; Publishing Company, 2nd Edition, p. 510
  26
  27 (define (newton:find-integer-root f df/dx x_0)
  28   (let loop ((x x_0) (fx (f x_0)))
  29     (cond
  30      ((zero? fx) x)
  31      (else
  32       (let ((df (df/dx x)))
  33         (cond
  34          ((zero? df) #f)                ; stuck at local min/max
  35          (else
  36           (let* ((delta (quotient (+ fx (quotient df 2)) df))
  37                  (next-x (cond ((not (zero? delta)) (- x delta))
  38                                ((positive? fx) (- x 1))
  39                                (else (- x -1))))
  40                  (next-fx (f next-x)))
  41             (cond ((>= (abs next-fx) (abs fx)) x)
  42                   (else (loop next-x next-fx)))))))))))
  43
  44 (define (integer-sqrt y)
  45   (newton:find-integer-root (lambda (x) (- (* x x) y))
  46                             (lambda (x) (* 2 x))
  47                             (ash 1 (quotient (integer-length y) 2))))
  48
  49 (define (newton:find-root f df/dx x_0 prec)
  50   (if (and (negative? prec) (integer? prec))
  51       (let loop ((x x_0) (fx (f x_0)) (count prec))
  52         (cond ((zero? count) x)
  53               (else (let ((df (df/dx x)))
  54                       (cond ((zero? df) #f) ; stuck at local min/max
  55                             (else (let* ((next-x (- x (/ fx df)))
  56                                          (next-fx (f next-x)))
  57                                     (cond ((= next-x x) x)
  58                                           ((> (abs next-fx) (abs fx)) #f)
  59                                           (else (loop next-x next-fx
  60                                                       (+ 1 count)))))))))))
  61       (let loop ((x x_0) (fx (f x_0)))
  62         (cond ((< (abs fx) prec) x)
  63               (else (let ((df (df/dx x)))
  64                       (cond ((zero? df) #f) ; stuck at local min/max
  65                             (else (let* ((next-x (- x (/ fx df)))
  66                                          (next-fx (f next-x)))
  67                                     (cond ((= next-x x) x)
  68                                           ((> (abs next-fx) (abs fx)) #f)
  69                                           (else (loop next-x next-fx))))))))))))
  70
  71 ;;; H. J. Orchard, "The Laguerre Method for Finding the Zeros of
  72 ;;; Polynomials", IEEE Transactions on Circuits and Systems, Vol. 36,
  73 ;;; No. 11, November 1989, pp 1377-1381.
  74
  75 (define (laguerre:find-root f df/dz ddf/dz^2 z_0 prec)
  76   (if (and (negative? prec) (integer? prec))
  77       (let loop ((z z_0) (fz (f z_0)) (count prec))
  78         (cond ((zero? count) z)
  79               (else
  80                (let* ((df (df/dz z))
  81                       (ddf (ddf/dz^2 z))
  82                       (disc (sqrt (- (* df df) (* fz ddf)))))
  83                  (if (zero? disc)
  84                      #f
  85                      (let* ((next-z
  86                              (- z (/ fz (if (negative? (+ (* (real-part df)
  87                                                              (real-part disc))
  88                                                           (* (imag-part df)
  89                                                              (imag-part disc))))
  90                                             (- disc) disc))))
  91                             (next-fz (f next-z)))
  92                        (cond ((>= (magnitude next-fz) (magnitude fz)) z)
  93                              (else (loop next-z next-fz (+ 1 count))))))))))
  94       (let loop ((z z_0) (fz (f z_0)) (delta-z #f))
  95         (cond ((< (magnitude fz) prec) z)
  96               (else
  97                (let* ((df (df/dz z))
  98                       (ddf (ddf/dz^2 z))
  99                       (disc (sqrt (- (* df df) (* fz ddf)))))
 100                  ;;(print 'disc disc)
 101                  (if (zero? disc)
 102                      #f
 103                      (let* ((next-z
 104                              (- z (/ fz (if (negative? (+ (* (real-part df)
 105                                                              (real-part disc))
 106                                                           (* (imag-part df)
 107                                                              (imag-part disc))))
 108                                             (- disc) disc))))
 109                             (next-delta-z (magnitude (- next-z z))))
 110                        ;;(print 'next-z next-z )
 111                        ;;(print '(f next-z) (f next-z))
 112                        ;;(print 'delta-z delta-z 'next-delta-z next-delta-z)
 113                        (cond ((zero? next-delta-z) z)
 114                              ((and delta-z (>= next-delta-z delta-z)) z)
 115                              (else
 116                               (loop next-z (f next-z) next-delta-z)))))))))))
 117
 118 (define (laguerre:find-polynomial-root deg f df/dz ddf/dz^2 z_0 prec)
 119   (if (and (negative? prec) (integer? prec))
 120       (let loop ((z z_0) (fz (f z_0)) (count prec))
 121         (cond ((zero? count) z)
 122               (else
 123                (let* ((df (df/dz z))
 124                       (ddf (ddf/dz^2 z))
 125                       (tmp (* (+ deg -1) df))
 126                       (sqrt-H (sqrt (- (* tmp tmp) (* deg (+ deg -1) fz ddf))))
 127                       (df+sqrt-H (+ df sqrt-H))
 128                       (df-sqrt-H (- df sqrt-H))
 129                       (next-z
 130                        (- z (/ (* deg fz)
 131                                (if (>= (magnitude df+sqrt-H)
 132                                        (magnitude df-sqrt-H))
 133                                    df+sqrt-H
 134                                    df-sqrt-H)))))
 135                  (loop next-z (f next-z) (+ 1 count))))))
 136       (let loop ((z z_0) (fz (f z_0)))
 137         (cond ((< (magnitude fz) prec) z)
 138               (else
 139                (let* ((df (df/dz z))
 140                       (ddf (ddf/dz^2 z))
 141                       (tmp (* (+ deg -1) df))
 142                       (sqrt-H (sqrt (- (* tmp tmp) (* deg (+ deg -1) fz ddf))))
 143                       (df+sqrt-H (+ df sqrt-H))
 144                       (df-sqrt-H (- df sqrt-H))
 145                       (next-z
 146                        (- z (/ (* deg fz)
 147                                (if (>= (magnitude df+sqrt-H)
 148                                        (magnitude df-sqrt-H))
 149                                    df+sqrt-H
 150                                    df-sqrt-H)))))
 151                  (loop next-z (f next-z))))))))
 152
 153 (define (secant:find-root-1 f x0 x1 prec must-bracket?)
 154   (letrec ((stop?
 155             (cond ((procedure? prec) prec)
 156                   ((and (integer? prec) (negative? prec))
 157                    (lambda (x0 x1 fmax count)
 158                      (>= count (- prec))))
 159                   (else
 160                    (lambda (x0 f0 x1 f1 count)
 161                      (and (< (abs f0) prec)
 162                           (< (abs f1) prec))))))
 163            (bracket-iter
 164             (lambda (xlo flo glo xhi fhi ghi count)
 165               (define (step xnew fnew)
 166                 (cond ((or (= xnew xlo)
 167                            (= xnew xhi))
 168                        (let ((xmid (+ xlo (* 1/2 (- xhi xlo)))))
 169                          (if (= xnew xmid)
 170                              xmid
 171                              (step xmid (f xmid)))))
 172                       ((positive? fnew)
 173                        (bracket-iter xlo flo (if glo (* 0.5 glo) 1)
 174                                      xnew fnew #f
 175                                      (+ count 1)))
 176                       (else
 177                        (bracket-iter xnew fnew #f
 178                                      xhi fhi (if ghi (* 0.5 ghi) 1)
 179                                      (+ count 1)))))
 180               (if (stop? xlo flo xhi fhi count)
 181                   (if (> (abs flo) (abs fhi)) xhi xlo)
 182                   (let* ((fflo (if glo (* glo flo) flo))
 183                          (ffhi (if ghi (* ghi fhi) fhi))
 184                          (del (- (/ fflo (- ffhi fflo))))
 185                          (xnew (+ xlo (* del (- xhi xlo))))
 186                          (fnew (f xnew)))
 187                     (step xnew fnew))))))
 188     (let ((f0 (f x0))
 189           (f1 (f x1)))
 190       (cond ((<= f0 0 f1)
 191              (bracket-iter x0 f0 #f x1 f1 #f 0))
 192             ((<= f1 0 f0)
 193              (bracket-iter x1 f1 #f x0 f0 #f 0))
 194             (must-bracket? #f)
 195             (else
 196              (let secant-iter ((x0 x0)
 197                                (f0 f0)
 198                                (x1 x1)
 199                                (f1 f1)
 200                                (count 0))
 201                (cond ((stop? x0 f0 x1 f1 count)
 202                       (if (> (abs f0) (abs f1)) x1 x0))
 203                      ((<= f0 0 f1)
 204                       (bracket-iter x0 f0 #f x1 f1 #f count))
 205                      ((>= f0 0 f1)
 206                       (bracket-iter x1 f1 #f x0 f0 #f count))
 207                      ((= f0 f1) #f)
 208                      (else
 209                       (let* ((xnew (+ x0 (* (- (/ f0 (- f1 f0))) (- x1 x0))))
 210                              (fnew (f xnew))
 211                              (fmax (max (abs f1) (abs fnew))))
 212                         (secant-iter x1 f1 xnew fnew (+ count 1)))))))))))
 213
 214 (define (secant:find-root f x0 x1 prec)
 215   (secant:find-root-1 f x0 x1 prec #f))
 216 (define (secant:find-bracketed-root f x0 x1 prec)
 217   (secant:find-root-1 f x0 x1 prec #t))