#lang racket/base
(require racket/contract)
(provide/contract 
 (chebyshev-integral (-> natural-number/c natural-number/c procedure?))
 (eval-over-range (-> number? number? procedure? number?))
 (x^y (-> number? number? number?)))

; f(x,p,q) = int x^p (1-x)^q dx
(define (chebyshev-integral p q)
  (cond 
    [(= p 0) ; f = int (1-x)^q dx
     (let [(q+1 (add1 q))]
       (lambda (x)
         (- (/ (x^y (- 1 x) q+1)
               q+1))))]
    [(= q 0) ; f = int x^p dx
     (let [(p+1 (add1 p))]
       (lambda (x)
         (/ (x^y x p+1) p+1)))]
    [else 
     ; we do integration by parts, int u dv = uv - int v du
     ; let u =   (1-x)^q       dv = x^p dx
     ;    du = -q(1-x)^(q-1)dx  v = x^(p+1)/(p+1)
     ; this results in another integral of the same form,
     ; allowing an interesting recursion
     (let [(p+1 (add1 p))]
       (lambda (x)
         (/ (+ (* (x^y (- 1 x) q)
                  (x^y x p+1))
               (* q
                  ; The tricky part - we incorporate the lambda
                  ; returned by the function... into the lambda
                  ; returned by the function!
                  (apply (chebyshev-integral p+1 (sub1 q))
                         (list x))))
            p+1)))]))

(define (eval-over-range x y f)
  (- (f y) (f x)))

(define (x^y x y)
  (if (= x 0)
      (if (> y 0) 0
          (error (format "(x^y ~a ~a)") x y))
      (exp (* y (log x)))))