(module math.polynomial.chebyshev (chebyshev-poly
                                   chebyshev-poly?
                                   chebyshev-poly-min
                                   chebyshev-poly-max
                                   chebyshev-poly-coefs

                                   chebyshev-poly-order
                                   chebyshev-poly-convert

                                   build-chebyshev-poly
                                   build-chebyshev-fppoly
                                   ;; build-chebyshev-bfpoly

                                   chebyshev-poly-fun
                                   chebyshev-fppoly-fun
                                   ;; chebyshev-bfpoly-fun

                                   inline-chebyshev-fppoly-fun
                                   chebyshev-iter)
  (import scheme
          chicken.base
          chicken.type
          chicken.flonum
          srfi-133
          miscmacros
          math.base
          ;; math.racket.bigfloat
          math.flonum.functions
          math.racket-shim)

  (define-type chebyshev-poly (struct chebyshev-poly))

  (define-record chebyshev-poly
    min
    max
    coefs)

  (define chebyshev-poly make-chebyshev-poly)

  (: chebyshev-poly-order (chebyshev-poly -> number))
  (define (chebyshev-poly-order p)
    (define n (vector-length (chebyshev-poly-coefs p)))
    (if (zero? n) 0 (- n 1)))


  (: chebyshev-poly-convert (chebyshev-poly (any -> any) -> chebyshev-poly))
  (define (chebyshev-poly-convert p conv)
    (define mn (chebyshev-poly-min p))
    (define mx (chebyshev-poly-max p))
    (define cs (chebyshev-poly-coefs p))
    (chebyshev-poly (conv mn) (conv mx) (vector-map conv cs)))

  ;; (: make-build-chebyshev-poly
  ;;    (All (A) ((A A -> A) (A A -> A) (A A -> A) (A A -> A) (A -> A) (Integer -> A) (-> A)
  ;;              -> (A A Integer (A -> A) -> (chebyshev-poly A)))))
  (define-syntax-rule (make-build-chebyshev-poly num+ num- num* num/ numcos int->num pi-thnk)
    (lambda (mn mx n f)
      (define 0.num (int->num 0))
      (define 1.num (int->num 1))
      (define 2.num (int->num 2))
      (define 1/2.num (num/ 1.num 2.num))
      (define pi.num (pi-thnk))
      (define mx-mn (num- mx mn))
      (define mn+mx (num+ mn mx))
      (let ([n.num  (int->num n)])
        (define norm (num/ 2.num n.num))
        (define cs
          (let ((vec (make-vector n)))
            (let loop ((j 0))
              (unless (= j n)
                (vector-set! vec j
                             (let ([j  (int->num j)])
                               (num*
                                (let loop ((cj 0.num)
                                           (k 0))
                                  (if (= k n)
                                      cj
                                      (let* ((theta (num/ (num* pi.num (num+ (int->num k) 1/2.num)) n.num))
                                             (x (num* 1/2.num (num+ (num* mx-mn (numcos theta)) mn+mx))))
                                        (loop (add1 k)
                                              (num+ cj (num* (f x) (numcos (num* theta j))))))))
                                norm)))
                (loop (add1 j))))))
        (chebyshev-poly mn mx cs))))

  (define-syntax-rule (make-chebyshev-poly-fun num+ num- num* num/ int->num)
    (lambda (p)
      (define cs (chebyshev-poly-coefs p))
      (define n (vector-length cs))
      (define 0.num (int->num 0))
      (define 2.num (int->num 2))
      (cond
       [(zero? n)  (lambda (x) 0.num)]
       [else
        (let ((mn (chebyshev-poly-min p))
              (mx (chebyshev-poly-max p)))
          (lambda (x)
            (define i (- n 1))
            (define c (vector-ref cs i))
            (define y (num/ (num- (num* x 2.num) (num+ mn mx))
                            (num- mx mn)))
            (define y2 (num* y 2.num))
            (let loop ([i i]
                       [c c]
                       [d 0.num]
                       [dd 0.num])
              (cond [(zero? i)  (num+ (num* y d) (num- (num/ c 2.num) dd))]
                    [else
                     (let ([d   (num+ (num* y2 d) (num- c dd))]
                           [dd  d])
                       (loop (- i 1) (vector-ref cs (- i 1)) d dd))]))))])))

  (define build-chebyshev-poly
    (make-build-chebyshev-poly + - * / cos (lambda (x) x) (lambda () pi)))

  (define build-chebyshev-fppoly
    (make-build-chebyshev-poly fp+ fp- fp* fp/ fpcos fp (lambda () pi)))

  ;; (define build-chebyshev-bfpoly
  ;;   (make-build-chebyshev-poly Bigfloat bf+ bf- bf* bf/ bfcos bf (lambda () pi.bf)))

  (define chebyshev-poly-fun
    (make-chebyshev-poly-fun + - * / (lambda (x) x)))

  (define chebyshev-fppoly-fun
    (make-chebyshev-poly-fun fp+ fp- fp* fp/ fp))

  ;; (define chebyshev-bfpoly-fun
  ;;   (make-chebyshev-poly-fun Bigfloat bf+ bf- bf* bf/ bf))

  (define-syntax chebyshev-iter
    (syntax-rules ()
      [(_ y y2 d dd ())  d]
      [(_ y y2 d dd (c0))
       (fp+ (fp* y d) (fp- (fp/ c0 2.0) dd))]
      [(_ y y2 d dd (c0 c ...))
       (let ([d   (fp+ (fp* y2 d) (fp- c0 dd))]
             [dd  d])
         (chebyshev-iter y y2 d dd (c ...)))]))

  (define-syntax inline-chebyshev-fppoly-fun
    (er-macro-transformer
     (lambda (x r c)
       (let ((lower (cadr x))
             (upper (caddr x))
             (lst (reverse (cadddr x))))
         `(lambda (z)
            (define y (fp/ (fp- (fp+ z z) (fp+ ,lower ,upper))
                           (fp- ,upper ,lower)))
            (define y2 (fp+ y y))
            (let ([d   0.0]
                  [dd  0.0])
              (chebyshev-iter y y2 d dd ,(cdr lst)))))))))