#lang scheme/base
;;; PLT Scheme Science Collection
;;; chebyshev.ss
;;; Copyright (c) 2004-2008 M. Douglas Williams
;;;
;;; This library is free software; you can redistribute it and/or
;;; modify it under the terms of the GNU Lesser General Public
;;; License as published by the Free Software Foundation; either
;;; version 2.1 of the License, or (at your option) any later version.
;;;
;;; This library is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
;;; Lesser General Public License for more details.
;;;
;;; You should have received a copy of the GNU Lesser General Public
;;; License along with this library; if not, write to the Free
;;; Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
;;; 02111-1307 USA.
;;;
;;; -------------------------------------------------------------------
;;;
;;; This code in based on the Chebyshev Approximations in the GNU
;;; Scientific Library (GSL), which is licensed under the GPL.
;;;
;;; Version  Date      Description
;;; 0.1.0    08/07/04  This is the initial release of the Chebyshev
;;;                    approximations routines ported from GSL. (Doug
;;;                    Williams)
;;; 1.0.0    09/29/04  Added contracts for functions.  Marked as ready
;;;                    for Release 1.0.  (Doug Williams)
;;; 2.0.0    11/17/07  Added unchecked versions of functions and getting
;;;                    ready for PLT Scheme V4.0.  (Doug Williams)
;;; 3.0.0    06/09/08  Global changes for V4.0.  (Doug Williams)

(require "math.ss"
         (lib "contract.ss"))

;; Structure: chebyshev-series
(define-struct chebyshev-series
  (coefficients
   order
   (lower #:mutable)
   (upper #:mutable)))

;; Provide unchecked versions of procedures
(provide
 (rename-out (make-chebyshev-series-order unchecked-make-chebyshev-series-order)
             (chebyshev-series-init unchecked-chebyshev-series-init)
             (chebyshev-eval unchecked-chebyshev-eval)
             (chebyshev-eval-n unchecked-chebyshev-eval-n)))

;;; Define the contracts for the module.
(provide/contract
 (struct chebyshev-series ((coefficients (vectorof real?))
                           (order natural-number/c)
                           (lower real?)
                           (upper real?)))
 (make-chebyshev-series-order
  (-> natural-number/c chebyshev-series?))
 (chebyshev-series-init
  (-> chebyshev-series? procedure? real? real? void?))
 (chebyshev-eval
  (-> chebyshev-series? real? real?))
 (chebyshev-eval-n
  (-> chebyshev-series? natural-number/c real? real?)))

;;; make-chebyshev-series-order: real -> chebyshev-series
(define (make-chebyshev-series-order order)
  (make-chebyshev-series
   (make-vector (+ order 1))
   order
   0.0 0.0))

;;; chebyshev-series-init: chebyshev-series x function x real x real
;;;   -> void
(define (chebyshev-series-init cs func a b)
  (when (>= a b)
    (error 'chebyshev-series-init
           "null function interval [a,b]"))
  (set-chebyshev-series-lower! cs a)
  (set-chebyshev-series-upper! cs b)
  (let* ((order (chebyshev-series-order cs))
         (coefficients (chebyshev-series-coefficients cs))
         (bma (* 0.5 (- b a)))
         (bpa (* 0.5 (+ b a)))
         (fac (/ 2.0 (+ order 1.0)))
         (f (make-vector (+ order 1))))
    (do ((k 0 (+ k 1)))
        ((> k order) (void))
      (let ((y (cos (* pi (/ (+ k 0.5) (+ order 1))))))
        (vector-set! f k
                     (func (+ (* y bma) bpa)))))
    (do ((j 0 (+ j 1)))
        ((> j order) (void))
      (let ((sum 0.0))
        (do ((k 0 (+ k 1)))
          ((> k order) (void))
          (set! sum (+ sum (* (vector-ref f k)
                              (cos (/ (* pi j (+ k 0.5))
                                      (+ order 1)))))))
        (vector-set! coefficients j (* fac sum))))))

;;; chebyshev-eval: chebyshev-series x real -> real
;;; Evaluate the Chebyshev series at the given point.
(define (chebyshev-eval cs x)
  (let* ((coefs (chebyshev-series-coefficients cs))
         (order (chebyshev-series-order cs))
         (lower (chebyshev-series-lower cs))
         (upper (chebyshev-series-upper cs))
         (d 0.0)
         (dd 0.0)
         (y (/ (- (* 2.0 x) lower upper) (- upper lower)))
         (y2 (* 2.0 y)))
    (do ((j order (- j 1)))
        ((= j 0) (void))
      (let ((temp d))
        (set! d (+ (* y2 d) (- dd) (vector-ref coefs j)))
        (set! dd temp)))
    (+ (* y d) (- dd) (* 0.5 (vector-ref coefs 0)))))

;;; chebyshev-eval-n: chebyshev-series x integer x real -> real
;;; Evaluate the Chebyshev series at the given point.
(define (chebyshev-eval-n cs n x)
  (let* ((coefs (chebyshev-series-coefficients cs))
         (order (chebyshev-series-order cs))
         (lower (chebyshev-series-lower cs))
         (upper (chebyshev-series-upper cs))
         (d 0.0)
         (dd 0.0)
         (y (/ (- (* 2.0 x) lower upper) (- upper lower)))
         (y2 (* 2.0 y)))
    (do ((j (min n order) (- j 1)))
        ((= j 0) (void))
      (let ((temp d))
        (set! d (+ (* y2 d) (- dd) (vector-ref coefs j)))
        (set! dd temp)))
    (+ (* y d) (- dd) (* 0.5 (vector-ref coefs 0)))))