;last change 2009-03-16
;2009-03-02

(setq *read-default-float-format* 'long-float)

(defvar *quadrature-error* 1e-12)

(defun quadrature-simpson (f a b &optional (err *quadrature-error*))
  (labels ((simpson-rule
             (h y0 y1 y2)
             (* 1/3 h (+ y0 y2 (* 4 y1))))
           (quadrature-segment
             (x2 y0 y2 y4 h err sum-0-2-4)
             (let* ((x1 (- x2 h))
                    (x3 (+ x2 h))
                    (y1 (funcall f x1))
                    (y3 (funcall f x3))
                    (sum-0-1-2
                      (simpson-rule h y0 y1 y2))
                    (sum-2-3-4
                      (simpson-rule h y2 y3 y4))
                    (sum-0-1-2-3-4
                      (+ sum-0-1-2 sum-2-3-4)))
               (if (<= (abs (- sum-0-1-2-3-4 sum-0-2-4))
                       (* 15 err))
                 sum-0-1-2-3-4
                 (let ((h (/ h 2))
                       (err (/ err 2)))
                   (if (or (= h 0) (= err 0))
                     ;running into hardware limit
                     sum-0-1-2-3-4
                     (+ (quadrature-segment x1 y0 y1 y2
                                            h err
                                            sum-0-1-2)
                        (quadrature-segment x3 y2 y3 y4
                                            h err
                                            sum-2-3-4))))))))
    (let* ((h (/ (- b a) 2.0))
           (m (+ a h))
           (fa (funcall f a))
           (fb (funcall f b))
           (fm (funcall f m)))
      (quadrature-segment
        m
        fa fm fb
        (/ h 2)
        err
        (simpson-rule h fa fm fb)))))