; numerical integration

(define-syntax fold-of
  (syntax-rules (range in is)
    ((_ "z" f b e) (set! b (f b e)))
    ((_ "z" f b e (v range fst pst stp) c ...)
      (let* ((x fst) (p pst) (s stp)
             (le? (if (positive? s) <= >=)))
        (do ((v x (+ v s))) ((le? p v) b)
          (fold-of "z" f b e c ...))))
    ((_ "z" f b e (v range fst pst) c ...)
      (let* ((x fst) (p pst) (s (if (< x p) 1 -1)))
        (fold-of "z" f b e (v range x p s) c ...)))
    ((_ "z" f b e (v range pst) c ...)
      (fold-of "z" f b e (v range 0 pst) c ...))
    ((_ "z" f b e (x in xs) c ...)
      (do ((t xs (cdr t))) ((null? t) b)
        (let ((x (car t)))
          (fold-of "z" f b e c ...))))
    ((_ "z" f b e (x is y) c ...)
      (let ((x y)) (fold-of "z" f b e c ...)))
    ((_ "z" f b e p? c ...)
      (if p? (fold-of "z" f b e c ...)))
    ((_ f i e c ...)
      (let ((b i)) (fold-of "z" f b e c ...)))))

(define-syntax sum-of (syntax-rules ()
  ((_ arg ...) (fold-of + 0 arg ...))))

(define (int-rect f a b n)
  (let* ((w (/ (- b a) n))
         (mid (lambda (i) (+ a (/ w -2) (* i w))))
         (g (lambda (i) (* w (f (mid i))))))
    (sum-of (g i) (i range 1 (+ n 1)))))

(define (int-trap f a b n)
  (let* ((w (/ (- b a) n))
         (lo (lambda (i) (+ a (* (- i 1) w))))
         (hi (lambda (i) (+ a (* i w))))
         (g (lambda (i) (* w 1/2 (+ (f (lo i)) (f (hi i)))))))
    (sum-of (g i) (i range 1 (+ n 1)))))

(define (int-simp f a b n)
  (let* ((w (/ (- b a) n))
         (lo (lambda (i) (+ a (* (- i 1) w))))
         (mid (lambda (i) (+ a (/ w -2) (* i w))))
         (hi (lambda (i) (+ a (* i w))))
         (g (lambda (i) (* w 1/6 (+ (f (lo i)) (* 4 (f (mid i))) (f (hi i)))))))
    (sum-of (g i) (i range 1 (+ n 1)))))

(define (integrate method f a b . epsilon)
  (let ((g5 (method f a b 5))
        (g10 (method f a b 10))
        (mid (/ (+ a b) 2))
        (eps (if (null? epsilon) #e1e-7 (car epsilon))))
    (if (< (abs (- g10 g5)) eps) g10
      (+ (integrate method f a mid eps)
         (integrate method f mid b eps)))))

(define (cube x) (* x x x))
(display (int-rect cube 0 1 10000.)) (newline)
(display (int-trap cube 0 1 10000.)) (newline)
(display (int-simp cube 0 1 10000.)) (newline)

(define (approx-pi n)
  (inexact->exact
    (round
      (integrate int-simp
        (lambda (x) (/ (log x)))
          2 n 1e-7))))

(display (approx-pi 1e21)) (newline)

0.24999999875000073
0.2500000025
0.25
21127269486616129536