From bd9733926076885e3417b74de76e4c9c7bc56254 Mon Sep 17 00:00:00 2001 From: Bryan Newbold Date: Mon, 20 Feb 2017 00:05:28 -0800 Subject: Import Upstream version 2c7 --- root.scm | 66 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) (limited to 'root.scm') 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)) -- cgit v1.2.3