summaryrefslogtreecommitdiffstats
path: root/root.scm
diff options
context:
space:
mode:
authorBryan Newbold <bnewbold@robocracy.org>2017-02-20 00:05:28 -0800
committerBryan Newbold <bnewbold@robocracy.org>2017-02-20 00:05:28 -0800
commitbd9733926076885e3417b74de76e4c9c7bc56254 (patch)
tree2c99dced547d48407ad44cb0e45e31bb4d02ce43 /root.scm
parentfa3f23105ddcf07c5900de47f19af43d1db1b597 (diff)
downloadslib-bd9733926076885e3417b74de76e4c9c7bc56254.tar.gz
slib-bd9733926076885e3417b74de76e4c9c7bc56254.zip
Import Upstream version 2c7upstream/2c7
Diffstat (limited to 'root.scm')
-rw-r--r--root.scm66
1 files changed, 66 insertions, 0 deletions
diff --git a/root.scm b/root.scm
index 3c764a6..d561af6 100644
--- a/root.scm
+++ b/root.scm
@@ -149,3 +149,69 @@
df+sqrt-H
df-sqrt-H)))))
(loop next-z (f next-z))))))))
+
+(define (secant:find-root-1 f x0 x1 prec must-bracket?)
+ (letrec ((stop?
+ (cond ((procedure? prec) prec)
+ ((and (integer? prec) (negative? prec))
+ (lambda (x0 x1 fmax count)
+ (>= count (- prec))))
+ (else
+ (lambda (x0 f0 x1 f1 count)
+ (and (< (abs f0) prec)
+ (< (abs f1) prec))))))
+ (bracket-iter
+ (lambda (xlo flo glo xhi fhi ghi count)
+ (define (step xnew fnew)
+ (cond ((or (= xnew xlo)
+ (= xnew xhi))
+ (let ((xmid (+ xlo (* 1/2 (- xhi xlo)))))
+ (if (= xnew xmid)
+ xmid
+ (step xmid (f xmid)))))
+ ((positive? fnew)
+ (bracket-iter xlo flo (if glo (* 0.5 glo) 1)
+ xnew fnew #f
+ (+ count 1)))
+ (else
+ (bracket-iter xnew fnew #f
+ xhi fhi (if ghi (* 0.5 ghi) 1)
+ (+ count 1)))))
+ (if (stop? xlo flo xhi fhi count)
+ (if (> (abs flo) (abs fhi)) xhi xlo)
+ (let* ((fflo (if glo (* glo flo) flo))
+ (ffhi (if ghi (* ghi fhi) fhi))
+ (del (- (/ fflo (- ffhi fflo))))
+ (xnew (+ xlo (* del (- xhi xlo))))
+ (fnew (f xnew)))
+ (step xnew fnew))))))
+ (let ((f0 (f x0))
+ (f1 (f x1)))
+ (cond ((<= f0 0 f1)
+ (bracket-iter x0 f0 #f x1 f1 #f 0))
+ ((<= f1 0 f0)
+ (bracket-iter x1 f1 #f x0 f0 #f 0))
+ (must-bracket? #f)
+ (else
+ (let secant-iter ((x0 x0)
+ (f0 f0)
+ (x1 x1)
+ (f1 f1)
+ (count 0))
+ (cond ((stop? x0 f0 x1 f1 count)
+ (if (> (abs f0) (abs f1)) x1 x0))
+ ((<= f0 0 f1)
+ (bracket-iter x0 f0 #f x1 f1 #f count))
+ ((>= f0 0 f1)
+ (bracket-iter x1 f1 #f x0 f0 #f count))
+ ((= f0 f1) #f)
+ (else
+ (let* ((xnew (+ x0 (* (- (/ f0 (- f1 f0))) (- x1 x0))))
+ (fnew (f xnew))
+ (fmax (max (abs f1) (abs fnew))))
+ (secant-iter x1 f1 xnew fnew (+ count 1)))))))))))
+
+(define (secant:find-root f x0 x1 prec)
+ (secant:find-root-1 f x0 x1 prec #f))
+(define (secant:find-bracketed-root f x0 x1 prec)
+ (secant:find-root-1 f x0 x1 prec #t))