nintegrate   xtlang

(bind-func nintegrate
  "integrate using Simpsons rule - N must be odd"
  (lambda (f:[double,double]* xmin xmax N)
    (let ((sumeven:double 0.0)
          (sumodd:double 0.0)
          (x 0.0) (n 0)
          (interval (/ (- xmax xmin) (i64tod (- N 1)))))
      (set! n 2)
      (while (< n N)
        (set! x (+ xmin (* interval (i64tod (- n 1)))))
        (set! sumodd (+ sumodd (f x)))
        (set! n (+ n 2)))
      (set! n 3)
      (while (< n N)
        (set! x (+ xmin (* interval (i64tod (- n 1)))))
        (set! sumeven (+ sumeven (f x)))
        (set! n (+ n 2)))
      (* (+ (* sumeven 2.0) (* sumodd 4.0)
            (f xmin) (f xmax))
         (/ interval 3.0)))))