Import Upstream version 2c7upstream/2c7

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))