special-functions/exponential-integral.ss
#lang scheme/base
;;; PLT Scheme Science Collection
;;; special-functions/exponential-integral.ss
;;; Copyright (c) 2006-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 file implements the gamma functions.  It is based on the
;;; Special Functions in the GNU Scientific Library (GSL).
;;;
;;; Version  Date      Description
;;; 1.0.0    02/09/06  Initial implementation.  (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 (expint-E1 unchecked-expint-E1)
             (expint-E1-scaled unchecked-expint-E1-scaled)
             (expint-E2 unchecked-expint-E2)
             (expint-E2-scaled unchecked-expint-E2-scaled)
             (expint-Ei unchecked-expint-Ei)
             (expint-Ei-scaled unchecked-expint-Ei-scaled)))

(provide/contract
 (expint-E1
  (-> real? real?))
 (expint-E1-scaled
  (-> real? real?))
 (expint-E2
  (-> real? real?))
 (expint-E2-scaled
  (-> real? real?))
 (expint-Ei
  (-> real? real?))
 (expint-Ei-scaled
  (-> real? real?)))

(require "../machine.ss")
(require "../math.ss")
(require "../chebyshev.ss")

;;; Chebyshev series data
(define AE11-data
  '#( 0.121503239716065790
     -0.065088778513550150
      0.004897651357459670
     -0.000649237843027216
      0.000093840434587471
      0.000000420236380882
     -0.000008113374735904
      0.000002804247688663
      0.000000056487164441
     -0.000000344809174450
      0.000000058209273578
      0.000000038711426349
     -0.000000012453235014
     -0.000000005118504888
      0.000000002148771527
      0.000000000868459898
     -0.000000000343650105
     -0.000000000179796603
      0.000000000047442060
      0.000000000040423282
     -0.000000000003543928
     -0.000000000008853444
     -0.000000000000960151
      0.000000000001692921
      0.000000000000607990
     -0.000000000000224338
     -0.000000000000200327
     -0.000000000000006246
      0.000000000000045571
      0.000000000000016383
     -0.000000000000005561
     -0.000000000000006074
     -0.000000000000000862
      0.000000000000001223
      0.000000000000000716
     -0.000000000000000024
     -0.000000000000000201
     -0.000000000000000082
      0.000000000000000017))

(define AE11-cs
  (make-chebyshev-series
   AE11-data
   38 -1.0 1.0))

(define AE12-data
  '#( 0.582417495134726740
     -0.158348850905782750
     -0.006764275590323141
      0.005125843950185725
      0.000435232492169391
     -0.000143613366305483
     -0.000041801320556301
     -0.000002713395758640
      0.000001151381913647
      0.000000420650022012
      0.000000066581901391
      0.000000000662143777
     -0.000000002844104870
     -0.000000000940724197
     -0.000000000177476602
     -0.000000000015830222
      0.000000000002905732
      0.000000000001769356
      0.000000000000492735
      0.000000000000093709
      0.000000000000010707
     -0.000000000000000537
     -0.000000000000000716
     -0.000000000000000244
     -0.000000000000000058))

(define AE12-cs
  (make-chebyshev-series
   AE12-data
   24 -1.0 1.0))

(define E11-data
  '#(-16.11346165557149402600
       7.79407277874268027690
      -1.95540581886314195070
       0.37337293866277945612
      -0.05692503191092901938
       0.00721107776966009185
      -0.00078104901449841593
       0.00007388093356262168
      -0.00000620286187580820
       0.00000046816002303176
      -0.00000003209288853329
       0.00000000201519974874
      -0.00000000011673686816
       0.00000000000627627066
      -0.00000000000031481541
       0.00000000000001479904
      -0.00000000000000065457
       0.00000000000000002733
      -0.00000000000000000108))

(define E11-cs
  (make-chebyshev-series
   E11-data
   18 -1.0 1.0))

(define E12-data
  '#(-0.03739021479220279500
      0.04272398606220957700
     -0.13031820798497005440
      0.01441912402469889073
     -0.00134617078051068022
      0.00010731029253063780
     -0.00000742999951611943
      0.00000045377325690753
     -0.00000002476417211390
      0.00000000122076581374
     -0.00000000005485141480
      0.00000000000226362142
     -0.00000000000008635897
      0.00000000000000306291
     -0.00000000000000010148
      0.00000000000000000315))

(define E12-cs
  (make-chebyshev-series
   E12-data
   15 -1.0 1.0))

(define AE13-data
  '#(-0.605773246640603460
     -0.112535243483660900
      0.013432266247902779
     -0.001926845187381145
      0.000309118337720603
     -0.000053564132129618
      0.000009827812880247
     -0.000001885368984916
      0.000000374943193568
     -0.000000076823455870
      0.000000016143270567
     -0.000000003466802211
      0.000000000758754209
     -0.000000000168864333
      0.000000000038145706
     -0.000000000008733026
      0.000000000002023672
     -0.000000000000474132
      0.000000000000112211
     -0.000000000000026804
      0.000000000000006457
     -0.000000000000001568
      0.000000000000000383
     -0.000000000000000094
      0.000000000000000023))

(define AE13-cs
  (make-chebyshev-series
   AE13-data
   24 -1.0 1.0))

(define AE14-data
  '#(-0.18929180007530170
     -0.08648117855259871
      0.00722410154374659
     -0.00080975594575573
      0.00010999134432661
     -0.00001717332998937
      0.00000298562751447
     -0.00000056596491457
      0.00000011526808397
     -0.00000002495030440
      0.00000000569232420
     -0.00000000135995766
      0.00000000033846628
     -0.00000000008737853
      0.00000000002331588
     -0.00000000000641148
      0.00000000000181224
     -0.00000000000052538
      0.00000000000015592
     -0.00000000000004729
      0.00000000000001463
     -0.00000000000000461
      0.00000000000000148
     -0.00000000000000048
      0.00000000000000016
     -0.00000000000000005))

(define AE14-cs
  (make-chebyshev-series
   AE14-data
   25 -1.0 1.0))

(define xmaxt (- log-double-min))
(define xmax (- xmaxt (log xmaxt)))

;;; Implementation for E1, allowing for scaling by exp(x)
(define (expint-E1-impl x scale)
  (cond ((and (< x (- xmax))
              (not scale))
         +inf.0)
        ((<= x -10.0)
         (let ((s (* (/ 1.0 x) (if scale 1.0 (exp (- x)))))
               (r (unchecked-chebyshev-eval AE11-cs (+ (/ 20.0 x) 1.0))))
           (* s (+ 1.0 r))))
        ((<= x -4.0)
         (let ((s (* (/ 1.0 x) (if scale 1.0 (exp (- x)))))
               (r (unchecked-chebyshev-eval AE12-cs (/ (+ (/ 40.0 x) 7.0) 3.0))))
           (* s (+ 1.0 r))))
        ((<= x -1.0)
         (let ((ln-term (- (log (abs x))))
               (sf (if scale (exp x) 1.0))
               (r (unchecked-chebyshev-eval E11-cs (/ (+ (* 2.0 x) 5.0) 3.0))))
           (* sf (+ ln-term r))))
        ((= x 0.0)
         -inf.0)
        ((<= x 1.0)
         (let ((ln-term (- (log (abs x))))
               (sf (if scale (exp x) 1.0))
               (r (unchecked-chebyshev-eval E12-cs x)))
           (* sf (+ ln-term -0.6875 x r))))
        ((<= x 4.0)
         (let ((s (* (/ 1.0 x) (if scale 1.0 (exp (- x)))))
               (r (unchecked-chebyshev-eval AE13-cs (/ (- (/ 8.0 x) 5.0) 3.0))))
           (* s (+ 1.0 r))))
        ((or (<= x xmax)
             scale)
         (let ((s (* (/ 1.0 x) (if scale 1.0 (exp (- x)))))
               (r (unchecked-chebyshev-eval AE14-cs (- (/ 8.0 x) 1.0))))
           (* s (+ 1.0 r))))
        (else
         0.0)))

(define (expint-E2-impl x scale)
  (cond ((and (< x (- xmax))
              (not scale))
         +inf.0)
        ((< x 100.0)
         (let ((ex (if scale 1.0 (exp (- x))))
               (r (expint-E1-impl x scale)))
           (- ex (* x r))))
        ((or (< x xmax)
             scale)
         (let* ((s (if scale 1.0 (exp (- x))))
                (c1  -2.0)
                (c2   6.0)
                (c3  -24.0)
                (c4   120.0)
                (c5  -720.0)
                (c6   5040.0)
                (c7  -40320.0)
                (c8   362880.0)
                (c9  -3628800.0)
                (c10  39916800.0)
                (c11 -479001600.0)
                (c12  6227020800.0)
                (c13 -87178291200.0)
                (y (/ 1.0 x))
                (sum9 (+ c9
                         (* y
                            (+ c10
                               (* y
                                  (+ c11
                                     (* y
                                        (+ c12
                                           (* y c13)))))))))
                (sum5 (+ c5
                         (* y
                            (+ c6
                               (* y
                                  (+ c7
                                     (* y
                                        (+ c8
                                           (* y sum9)))))))))
                (sum (* y
                        (+ c1
                           (* y
                              (+ c2
                                 (* y
                                    (+ c3
                                       (* y
                                          (+ c4
                                             (* y sum5)))))))))))
           (* s (/ (+ 1.0 sum) x))))
        (else
         0.0)))

(define (expint-E1 x)
  (expint-E1-impl x #f))

(define (expint-E1-scaled x)
  (expint-E1-impl x #t))

(define (expint-E2 x)
  (expint-E2-impl x #f))

(define (expint-E2-scaled x)
  (expint-E2-impl x #t))

(define (expint-Ei x)
  (- (expint-E1 (- x))))

(define (expint-Ei-scaled x)
  (- (expint-E1-scaled (- x))))