对于lisp党来说, 一切控制结构都是纸老虎

;; -*- coding:utf-8; syntax:common-lisp -*-

(defun bisection-recur (f a b &optional (e 1.0L-6) (recur-nth 0))
  "Page 79. bisection"
  (let ((alpha (/ (+ a b) 2.0L0)))            ; (1) use long-float
    (format t "k=~D a=~F b=~F f(alpha)=~F |b-a|=~F~%"
            recur-nth a b (funcall f alpha) (abs (- b a)))
    (cond ((or (< (abs (- b alpha)) e) 
               (< (abs (funcall f alpha)) e))
            alpha)                      ; (2)
          ((< (* (funcall f alpha) 
                 (funcall f a)) 
              0)
           (bisection f a alpha e (1+ recur-nth))) ;(3)
          (t (bisection f alpha b e (1+ recur-nth)))))) ; (4)
          

(defun bisection (f a b &optional (e 1.0L-6)
                   &key (max-recur 100))
  "Page 79. bisection"
  (loop 
     :for k from 0 to max-recur
     :for alpha = (/ (+ a b) 2.0L0)     ; (1)
     :do (format t "k=~D a=~F b=~F f(alpha)=~F |b-a|=~F~%"
                 k a b (funcall f alpha) (abs (- b a)))
     :thereis (when (or (< (abs (- b alpha)) e) 
                        (< (abs (funcall f alpha)) e))
                alpha)                  ; (2)
     :when (< (* (funcall f alpha) (funcall f a)) 0) :do (setq b alpha) ; (3)
     :else :do (setq a alpha)))         ; (4)

;; TEST: (simple-iteration-aitken #'testphi 5)
(defun simple-iteration-aitken (phi x0 &optional (e 1.0L-6)
                                    &key (max-recur 100))
  "Page 87. Aitken speedup."
  (let ((x (coerce x0 'long-float))
        (alpha nil))
    (format t "x_~D=~F~%" 0 x)
    (loop
       :for k from 0 to max-recur
       :for y = (funcall phi x)
       :for z = (funcall phi y)
       :when (= 0 (+ (- z (* 2 y)) x)) :do (return x)
       :do (setq alpha (- x
                          (/ (expt (- y x) 2)
                             (+ (- z (* 2 y)) x))))
       :do (format t "x_~D=~F |x_~D-x_~D|=~F~%"
                   (1+ k) alpha (1+ k) k (abs (- alpha x)))
       :when (< (abs (- alpha x)) e) :do (return alpha)
       :do (setq x alpha))))
     

;; TEST: (newton-iteration #'testfun #'testfun^ 5)
(defun newton-iteration (f f^ x0 &optional (e 1.0L-6)
                           &key (max-recur 100))
  "Page 92. Newton iter."
  (let ((x (coerce x0 'long-float))
        (alpha nil))
    (loop 
       :for k from 0 to max-recur
       :for fx = (funcall f x)
       :for f^x = (funcall f^ x)        ; (2)
       :do (format t "k=~D x_~D=~F f(x_~D)=~F "
                   k k x k fx)
       :when (= f^x 0) :do (return nil) ; (3)
       :do (setq alpha (- x (/ fx f^x))) ; (4)
           (format t "|x_~D-x_~D|=~F ~%" 
                   (1+ k) k (abs (- alpha x)))
       :thereis (when (< (abs (- alpha x)) e) ; (5)
                  alpha)
       :do (setq x alpha))))            ; (7)


;; TEST: (newton-iteration-simplify #'sin 0.7 5)
(defun newton-iteration-simplify (f M x0 &optional (e 1.0L-6)
                                    &key (max-recur 100))
  "Page 92. Newton iter simplify."
  (let ((x (coerce x0 'long-float))
        (alpha nil))
    (loop 
       :for k from 0 to max-recur
       :for fx = (funcall f x)
       :do (format t "k=~D x_~D=~F f(x_~D)=~F "
                   k k x k fx)
           (setq alpha (- x (/ fx M)))
           (format t "|x_~D-x_~D|=~F ~%" 
                   (1+ k) k (abs (- alpha x)))
       :thereis (when (< (abs (- alpha x)) e)
                  alpha)
       :do (setq x alpha))))

;; TEST: (secant-method #'testfun 5 4)
(defun secant-method (f x0 x1 &optional (e 1.0L-6)
                        &key (max-recur 100)x)
  "Page 92. Newton iter simplify."
  (let ((xk_1 (coerce x0 'long-float))
        (xk (coerce x1 'long-float))
        (alpha nil))
    (loop 
       :for k from 1 to max-recur
       :for fx = (funcall f xk)
       :for f^x = (/ (- fx (funcall f xk_1))
                     (- xk xk_1))
       :do (format t "k=~D x_~D=~F f(x_~D)=~F "
                   k k xk k fx)
       :when (= f^x 0) :do (return nil)
       :do (setq alpha (- xk (/ fx f^x)))
           (format t "|x_~D-x_~D|=~F ~%" 
                   (1+ k) k (abs (- alpha xk)))
       :thereis (when (< (abs (- alpha xk)) e)
                  alpha)
       :do (setq xk_1 xk)
           (setq xk alpha))))




;; test case
        
(defun testfun (x)
  (- (* x (exp x)) 1))

(defun testfun^ (x)
  (+ (exp x)
     (* x (exp x))))

;; use alpha=pi/2 to test
(defun testphi (x)
  (+ 1.6 (* 0.99 (cos x))))