;;; PLT Scheme Science Collection
;;; special-functions/error.ss
;;; Copyright (c) 2004 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 Error Functions in the GNU Scientific
;;; Library (GSL).
;;;
;;; Version  Date      Description
;;; 0.1.0    08/05/04  This is the initial release of the guassion
;;;                    distribution routines ported from GSL. (Doug
;;;                    Williams)
;;; 1.0.0    09/28/04  Added log-erfc and hazard functions.  Marked as
;;                     ready for Release 1.0.  Added contracts for
;;;                    functions.  (Doug Williams)

(module error mzscheme
  
  (require (lib "contract.ss"))
  
  (provide/contract
   (erfc
    (-> real? (real-in 0.0 2.0)))
   (log-erfc
    (-> real? real?))
   (erf
    (-> real? (real-in -1.0 1.0)))
   (hazard
    (-> real? (>=/c 0.0))))
  
  (require "../machine.ss")
  (require "../math.ss")
  (require "../chebyshev.ss")

  ;; Estimates erfc(x) valid for 8 < x < 100.
  ;; This is based on index 5725 in Hart, et. al.

  (define (erfc8-sum x)
    (define P
      #(2.97886562639399288862
        7.409740605964741794425
        6.1602098531096305440906
        5.019049726784267463450058
        1.275366644729965952479585264
        0.5641895854377550741253201704))
    (define Q
      #(3.3690752069827527677
        9.608965327192787870698
        17.08144074746600431571095
        12.0489519278551290360340491
        9.396034016235054150430579648
        2.260528520767326969591866945
        1.0))
    (let ((num (vector-ref P 5))
          (den (vector-ref Q 6)))
      (do ((i 4 (- i 1)))
        ((< i 0) (void))
        (set! num (+ (* x num) (vector-ref P i))))
      (do ((i 5 (- i 1)))
        ((< i 0) (void))
        (set! den (+ (* x den) (vector-ref Q i))))
      (/ num den)))

  (define (erfc8 x)
    (let ((e (erfc8-sum x)))
      (* e (exp (* (- x) x)))))
  
  (define (log-erfc8 x)
    (let ((e (erfc8-sum x)))
      (- (log e) (* x x))))

  ;; Abramowitz+Stegun, 7.1.5
  
  (define (erfseries x)
    (let* ((coef x)
           (e coef)
           (del 0.0))
      (do ((k 1 (+ k 1)))
        ((= k 30) (void))
        (set! coef (* coef (/ (* (- x) x) k)))
        (set! del (/ coef (+ (* 2.0 k) 1.0)))
        (set! e (+ e del)))
      (* (/ 2.0 (sqrt pi)) e)))

  ;; Chebyshev fit for erfc((t+1)/2, -1 < t < 1
  
  (define erfc-xlt1-data
    #( 1.06073416421769980345174155056
      -0.42582445804381043569204735291
       0.04955262679620434040357683080
       0.00449293488768382749558001242
      -0.00129194104658496953494224761
      -0.00001836389292149396270416979
       0.00002211114704099526291538556
      -5.23337485234257134673693179020e-7      
      -2.78184788833537885382530989578e-7
       1.41158092748813114560316684249e-8
       2.72571296330561699984539141865e-9
      -2.06343904872070629406401492476e-10
      -2.14273991996785367924201401812e-11
       2.22990255539358204580285098119e-12
       1.36250074650698280575807934155e-13
      -1.95144010922293091898995913038e-14
      -6.85627169231704599442806370690e-16
       1.44506492869699938239521607493e-16
       2.45935306460536488037576200030e-18
      -9.29599561220523396007359328540e-19))
  
  (define erfc-xlt1-cs
    (make-chebyshev-series
     erfc-xlt1-data
     19 -1.0 1.0))

  ;; Chebyshev fit for erfc(x) exp(x^2), 1 < x < 5, x = 2t + 3,
  ;; -1 < t < 1
  
  (define erfc-x15-data
    #( 0.44045832024338111077637466616
      -0.143958836762168335790826895326
       0.044786499817939267247056666937
      -0.013343124200271211203618353102
       0.003824682739750469767692372556
      -0.001058699227195126547306482530
       0.000283859419210073742736310108
      -0.000073906170662206760483959432
       0.000018725312521489179015872934
      -4.62530981164919445131297264430e-6
       1.11558657244432857487884006422e-6
      -2.63098662650834130067808832725e-7
       6.07462122724551777372119408710e-8
      -1.37460865539865444777251011793e-8
       3.05157051905475145520096717210e-9
      -6.65174789720310713757307724790e-10
       1.42483346273207784489792999706e-10
      -3.00141127395323902092018744545e-11
       6.22171792645348091472914001250e-12
      -1.26994639225668496876152836555e-12
       2.55385883033257575402681845385e-13
      -5.06258237507038698392265499770e-14
       9.89705409478327321641264227110e-15
      -1.90685978789192181051961024995e-15
       3.50826648032737849245113757340e-16))

  (define erfc-x15-cs
    (make-chebyshev-series
     erfc-x15-data
     24 -1.0 1.0))
  
  ;; Chebyshev fit for erfc(x) x exp(x^2), 5 < x < 10,
  ;; x = (5t + 15)/2, -1 < t < 1
  
  (define erfc-x510-data
    #( 1.11684990123545698684297865808
       0.003736240359381998520654927536
      -0.000916623948045470238763619870
       0.000199094325044940833965078819
      -0.000040276384918650072591781859
       7.76515264697061049477127605790e-6
      -1.44464794206689070402099225301e-6
       2.61311930343463958393485241947e-7
      -4.61833026634844152345304095560e-8
       8.00253111512943601598732144340e-9
      -1.36291114862793031395712122089e-9
       2.28570483090160869607683087722e-10
      -3.78022521563251805044056974560e-11
       6.17253683874528285729910462130e-12
      -9.96019290955316888445830597430e-13
       1.58953143706980770269606726000e-13
      -2.51045971047162509999527428316e-14
       3.92607828989125810013581287560e-15
      -6.07970619384160374392535453420e-16
       9.12600607264794717315507477670e-17))
  
  (define erfc-x510-cs
    (make-chebyshev-series
     erfc-x510-data
     19 -1.0 1.0))

  ;; Complementary error function, erfc
  ;;
  ;; This function computes the complementary error function of x.
  (define (erfc x)
    (let* ((ax (abs x))
           (val 
            (cond ((<= ax 1.0)
                   (let ((t (- (* 2.0 ax) 1.0)))
                     (chebyshev-eval erfc-xlt1-cs t)))
                  ((<= ax 5.0)
                   (let ((ex2 (exp (* (- x) x)))
                         (t (* 0.5 (- ax 3.0))))
                     (* ex2 
                        (chebyshev-eval erfc-x15-cs t))))
                  ((<= ax 10.0)
                   (let ((exterm (/ (exp (* (- x) x)) ax))
                         (t (/ (- (* 2.0 ax) 15.0) 5.0)))
                     (* exterm 
                        (chebyshev-eval erfc-x510-cs t))))
                  (else
                   (erfc8 ax)))))
      (if (< x 0.0)
          (- 2.0 val)
          val)))
  
  ;; log_erfc: real -> real
  ;;
  ;; This function computes the logarithm of the complementary error
  ;; function.
  (define (log-erfc x)
    (cond ((< (* x x) root6-double-epsilon)
           (let* ((y (/ x sqrtpi))
                  ;; series for -1/2 log(erfc(sqrt(ppi) y))
                  (c3 (/ (- 4.0 pi) 3.0))
                  (c4 (* 2.0 (/ (- 1.0 pi) 3.0)))
                  (c5 -0.001829764677455021)
                  (c6  0.02629651521057465)
                  (c7 -0.01621575378835404)
                  (c8  0.00125993961762116)
                  (c9  0.00556964649138)
                  (c10 -0.0045563339802)
                  (c11  0.0009461589032)
                  (c12  0.0013200243174)
                  (c13 -0.00142906)
                  (c14  0.00048204)
                  (series (+ c8 (* y (+ c9 (* y (+ c10 (* y (+ c11
                            (* y (+ c12 (* y (+ c12 (* y c14))))))))))))))
             (set! series (* y (+ 1.0 (* y (+ 1.0 (* y (+ c3 (* y (+ c4
                            (* y (+ c5 (* y (+ c6 (* y (+ c7 (* y series))))))))))))))))
             (* -2.0 series)))
          ((< (abs x) 1.0)
           (log1p (- (erf x))))
          ((> x 8.0)
           (log-erfc8 x))
          (else
           (log (erfc x)))))

  ;; erf: real -> real
  ;; This function computes the error function.
  (define (erf x)
    (if (< (abs x) 1.0)
        (erfseries x)
        (- 1.0 (erfc x))))
  
  ;; hazard: real -> real
  ;; This function computes the hazard function for the normal
  ;; distribution.
  (define (hazard x)
    (if (< x 25.0)
        (let* ((ln-erfc-val (log-erfc (/ x sqrt2)))
               (lnc -0.22579135264472743236)
               (arg (- lnc (* 0.5 x x) ln-erfc-val)))
          (exp arg))
        (let* ((ix2 (/ 1.0 (* x x)))
               (corrB (- 1.0 (* 9.0 ix2 (- 1.0 (* 11.0 ix2)))))
               (corrM (- 1.0 (* 5.0 ix2) (- 1.0 (* 7.0 ix2 corrB))))
               (corrT (- 1.0 (* ix2 (- 1.0 (* 3.0 ix2 corrM))))))
          (/ x corrT))))
    
)