;Dorai Sitaram
;ds26gte at yahoo dot com
;2009-03-02
;last change 2009-03-16

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

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

(defun quadrature-boole (f a b &optional (err *quadrature-error*))
  (labels ((boole-rule
             (h y0 y1 y2 y3 y4)
             (* 2/45 h (+ (* 7 y0)
                          (* 32 y1)
                          (* 12 y2)
                          (* 32 y3)
                          (* 7 y4))))
           (quadrature-segment
             (x1 x3
                 y0 y1 y2 y3 y4
                 h err sum4)
             (let* ((x0.5 (- x1 h))
                    (x1.5 (+ x1 h))
                    (x2.5 (- x3 h))
                    (x3.5 (+ x3 h))
                    (y0.5 (funcall f x0.5))
                    (y1.5 (funcall f x1.5))
                    (y2.5 (funcall f x2.5))
                    (y3.5 (funcall f x3.5))
                    (sum4-left
                      (boole-rule h y0 y0.5 y1 y1.5 y2))
                    (sum4-right
                      (boole-rule h y2 y2.5 y3 y3.5 y4))
                    (sum4-left-plus-sum4-right
                      (+ sum4-left sum4-right)))
               (if (<= (abs (- sum4-left-plus-sum4-right sum4))
                       (* 63 err))
                 sum4-left-plus-sum4-right
                 (let ((h (/ h 2))
                       (err (/ err 2)))
                   (if (or (= h 0) (= err 0))
                     ;running into hardware limit
                     sum4-left-plus-sum4-right
                     (+ (quadrature-segment
                          x0.5  x1.5
                          y0 y0.5 y1 y1.5 y2
                          h err
                          sum4-left)
                        (quadrature-segment
                          x2.5 x3.5
                          y2 y2.5 y3 y3.5 y4
                          h err
                          sum4-right))))))))
    (let* ((h (/ (- b a) 4.0))
           (x1 (+ a h))
           (x3 (- b h))
           (y0 (funcall f a))
           (y1 (funcall f x1))
           (y2 (funcall f (+ x1 h)))
           (y3 (funcall f x3))
           (y4 (funcall f b)))
      (quadrature-segment
        x1 x3
        y0 y1 y2 y3 y4
        (/ h 2)
        err
        (boole-rule h y0 y1 y2 y3 y4)))))

(defun quadrature (f a b &optional (err *quadrature-error*))
  (labels ((quadrature-aux
             (a b err)
             (cond ((and (not a) (not b))
                    (let ((err (/ err 2)))
                      (+ (quadrature-aux nil 0 err)
                         (quadrature-aux 0 nil err))))
                   ((not a)
                    (let ((f2 (lambda (u)
                                (if (= u 0) 0
                                  (/ (funcall f (log u))
                                     u)))))
                      (quadrature-boole f2 0 (exp b) err)))
                   ((not b)
                    (let ((f2 (lambda (u)
                                (if (= u 0) 0
                                  (/ (funcall f (- (log u)))
                                     u)))))
                      (quadrature-boole f2 0 (exp (- a)) err))))))
    (if (and a b)
      (quadrature-boole f a b err)
      (quadrature-aux a b err))))