Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
branch: master
Fetching contributors…

Octocat-spinner-32-eaf2f5

Cannot retrieve contributors at this time

file 90 lines (86 sloc) 2.948 kb
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
(defun one-newton-step (x f df)
  (let* ((fx (pull-bits (funcall f x)))
         (dfx (pull-bits (funcall df x))))
    (values (and (/= dfx 0)
                 (pull-bits (- x (/ fx dfx))))
            fx x)))

;; try to bracket a root of f.
(defun initial-newton (lo hi f df)
  (declare (optimize debug))
  (let* ((*precision* 64)
         (lo (round-to-float lo))
         (hi (round-to-float hi))
         (slo (signum (pull-bits (funcall f lo))))
         (shi (signum (pull-bits (funcall f hi))))
         (x (/ (+ lo hi) 2)))
    (block nil
      (cond ((zerop slo)
             (return (values lo lo)))
            ((zerop shi)
             (return (values hi hi)))
            ((/= slo shi)
             (return (values lo hi))))
      (loop repeat 32 do
       (multiple-value-bind (x2 fx)
           (one-newton-step (pull-bits x) f df)
         (cond ((/= (signum fx) slo)
                (return
                  (if (and x2 (< x2 x))
                      (values lo x)
                      (values x hi))))
               ((not x2)
                (return))
               ((>= x2 hi)
                (return))
               ((<= x2 lo)
                (return))
               (t (setf x x2))))))))

(defun bracketed-newton (lo hi f df)
  (let* ((*precision* 64)
         (lo (pull-bits lo))
         (hi (pull-bits hi))
         (slo (signum (pull-bits (funcall f lo))))
         (shi (signum (pull-bits (funcall f hi))))
         (x (/ (+ lo hi) 2)))
    (block nil
      (cond ((zerop slo)
             (return lo))
            ((zerop shi)
             (return hi))
            ((<= (abs (- (float-bits lo)
                         (float-bits hi)))
                 1)
             (return x)))
      (loop for i upfrom 1 do
        (let ((*precision* (min (max (ash 1 i) 64)
                                1024))
              (newtonp t))
          (multiple-value-bind (x2 fx)
              (and (plusp (mod i 3))
                   (one-newton-step (pull-bits x) f df))
            (when (or (not x2)
                      (>= x2 (max hi lo))
                      (<= x2 (min hi lo)))
              (setf newtonp nil
                    x (/ (+ lo hi) 2)
                    x2 x
                    fx (pull-bits (funcall f x))))
            (when (zerop fx)
              (return x))
            (let ((sx (signum fx)))
              (cond ((= sx shi)
                     (setf hi x))
                    ((= sx slo)
                     (setf lo x))))
            (setf x (if newtonp
                        (pull-bits x2)
                        (/ (+ lo hi) 2)))
            (when (<= (abs (- (float-bits lo)
                              (float-bits hi)))
                      1)
              (return x))))))))

(defun newton (lo hi f df)
  (multiple-value-bind (lo hi)
      (initial-newton (round-to-float lo) (round-to-float hi) f df)
    (and lo hi
         (bracketed-newton lo hi f df))))
Something went wrong with that request. Please try again.