#lang scheme/base
;;; PLT Scheme Science Collection
;;; special-functions/error.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 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)
;;; 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  Changes required for V4.0.  (Doug Williams)

(require (lib "contract.ss"))

(provide
 (rename-out (erfc unchecked-erfc)
             (log-erfc unchecked-log-erfc)
             (erf unchecked-erf)
             (hazard unchecked-hazard)))

(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)))
                   (unchecked-chebyshev-eval erfc-xlt1-cs t)))
                ((<= ax 5.0)
                 (let ((ex2 (exp (* (- x) x)))
                       (t (* 0.5 (- ax 3.0))))
                   (* ex2 
                      (unchecked-chebyshev-eval erfc-x15-cs t))))
                ((<= ax 10.0)
                 (let ((exterm (/ (exp (* (- x) x)) ax))
                       (t (/ (- (* 2.0 ax) 15.0) 5.0)))
                   (* exterm 
                      (unchecked-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(pi) 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))))