1# exponential.tcl --
2#    Compute exponential integrals (E1, En, Ei, li, Shi, Chi, Si, Ci)
3#
4
5namespace eval ::math::special {
6    variable pi 3.1415926
7    variable gamma 0.57721566490153286
8    variable halfpi [expr {$pi/2.0}]
9
10# Euler's digamma function for small integer arguments
11
12    variable psi {
13        NaN
14        -0.57721566490153286 0.42278433509846713 0.92278433509846713
15        1.2561176684318005 1.5061176684318005 1.7061176684318005
16        1.8727843350984672 2.0156414779556102 2.1406414779556102
17        2.2517525890667214 2.3517525890667215 2.4426616799758123
18        2.5259950133091458 2.6029180902322229 2.6743466616607945
19        2.7410133283274614 2.8035133283274614 2.8623368577392259
20        2.9178924132947812 2.9705239922421498 3.0205239922421496
21        3.0681430398611971 3.1135975853157425 3.1570758461853079
22        3.1987425128519744 3.2387425128519745 3.2772040513135128
23        3.31424108835055 3.3499553740648356 3.3844381326855251
24        3.4177714660188583 3.4500295305349873 3.4812795305349873
25        3.5115825608380176 3.5409943255438998 3.5695657541153283
26        3.597343531893106 3.6243705589201332 3.6506863483938172
27        3.6763273740348428
28    }
29}
30
31# ComputeExponFG --
32#    Compute the auxiliary functions f and g
33#
34# Arguments:
35#    x            Parameter of the integral (x>=0)
36# Result:
37#    Approximate values for f and g
38# Note:
39#    See Abramowitz and Stegun
40#
41proc ::math::special::ComputeExponFG {x} {
42    set x2 [expr {$x*$x}]
43    set fx [expr {($x2*$x2+7.241163*$x2+2.463936)/
44                  ($x2*$x2+9.068580*$x2+7.157433)/$x}]
45    set gx [expr {($x2*$x2+7.547478*$x2+1.564072)/
46                  ($x2*$x2+12.723684*$x2+15.723606)/$x2}]
47    list $fx $gx
48}
49
50
51# exponential_Ei --
52#    Compute the exponential integral of the second kind, to relative
53#    error eps
54# Arguments:
55#    x       Value of the argument
56#    eps     Relative error
57# Result:
58#    Principal value of the integral exp(x)/x
59#    from -infinity to x
60#
61proc ::math::special::exponential_Ei { x { eps 1.0e-10 } } {
62    variable gamma
63
64    if { ![string is double -strict $x] } {
65        return -code error "expected a floating point number but found \"$x\""
66    }
67    if { $x < 0.0 } {
68        return [expr { -[exponential_En 1 [expr { - $x }] $eps] }]
69    }
70    if { $x == 0.0 } {
71       set message "Argument to exponential_Ei must not be zero"
72       return -code error -errorcode [list ARITH DOMAIN $message] $message
73    }
74    if { $x >= -log($eps) } {
75        # evaluate Ei(x) as an asymptotic series; the series is formally
76        # divergent, but the leading terms give the desired value to
77        # enough precision.
78        set sum 0.
79        set term 1.
80        set k 1
81        while { 1 } {
82            set p $term
83            set term [expr { $term * ( $k / $x ) }]
84            if { $term < $eps } {
85                break
86            }
87            if { $term < $p } {
88                set sum [expr { $sum + $term }]
89            } else {
90                set sum [expr { $sum - $p }]
91                break
92            }
93            incr k
94        }
95        return [expr { exp($x) * ( 1.0 + $sum ) / $x }]
96    } elseif { $x >= 1e-18 } {
97        # evaluate Ei(x) as a power series
98        set sum 0.
99        set fact 1.
100        set pow $x
101        set n 1
102        while { 1 } {
103            set fact [expr { $fact * $n }]
104            set term [expr { $pow / $n / $fact }]
105            set sum [expr { $sum + $term }]
106            if { $term < $eps * $sum } break
107            set pow [expr { $x * $pow }]
108            incr n
109        }
110        return [expr { $sum + $gamma + log($x) }]
111    } else {
112        # Ei(x) for small x
113        return [expr { log($x) + $gamma }]
114    }
115}
116
117
118# exponential_En --
119#    Compute the exponential integral of n-th order, to relative error
120#    epsilon
121#
122# Arguments:
123#    n            Order of the integral (n>=1, integer)
124#    x            Parameter of the integral (x>0)
125#    epsilon      Relative error
126# Result:
127#    Value of En(x) = integral from 0 to x of exp(-x)/x**n
128#
129proc ::math::special::exponential_En { n x { epsilon 1.0e-10 } } {
130    variable psi
131    variable gamma
132    if { ![string is integer -strict $n] || $n < 0 } {
133        return -code error "expected a non-negative integer but found \"$n\""
134    }
135    if { ![string is double -strict $x] } {
136        return -code error "expected a floating point number but found \"$x\""
137    }
138    if { $n == 0 } {
139        if { $x == 0.0 } {
140            return -code error "E0(0) is indeterminate"
141        }
142        return [expr { exp( -$x ) / $x }]
143    }
144    if { $n == 1 && $x < 0.0 } {
145        return [expr {- [exponential_Ei [expr { -$x }] $eps] }]
146    }
147    if { $x < 0.0 } {
148        return -code error "can't evaluate En(x) for negative x"
149    }
150    if { $x == 0.0 } {
151        return [expr { 1.0 / ( $n - 1 ) }]
152    }
153
154    if { $x > 1.0 } {
155        # evaluate En(x) as a continued fraction
156        set b [expr { $x + $n }]
157        set c 1.e308
158        set d [expr { 1.0 / $b }]
159        set h $d
160        set i 1
161        while { 1 } {
162            set a [expr { -$i * ( $n - 1 + $i ) }]
163            set b [expr { $b + 2.0 }]
164            set d [expr { 1.0 / ( $a * $d + $b ) }]
165            set c [expr { $b + $a / $c }]
166            set delta [expr { $c * $d }]
167            set h [expr { $h * $delta }]
168            if { abs( $delta - 1. ) < $epsilon } {
169                return [expr { $h * exp(-$x) }]
170            }
171            incr i
172        }
173    } else {
174        # evaluate En(x) as a series
175        if { $n == 1 } {
176            set a [expr { -log($x) - $gamma }]
177        } else {
178            set a [expr { 1.0 / ( $n - 1 ) }]
179        }
180        set f 1.0
181        set i 1
182        while { 1 } {
183            set f [expr { - $f *  $x / $i }]
184            if { $i == $n - 1 } {
185                set term [expr { $f * ([lindex $psi $n] - log($x)) }]
186            } else {
187                set term [expr { $f / ( $n - 1 - $i ) }]
188            }
189            set a [expr { $a + $term }]
190            if { abs($term) < $epsilon * abs($a) } {
191                return $a
192            }
193            incr i
194        }
195    }
196}
197
198# exponential_E1 --
199#    Compute the exponential integral
200#
201# Arguments:
202#    x            Parameter of the integral (x>0)
203# Result:
204#    Value of E1(x) = integral from x to infinity of exp(-x)/x
205# Note:
206#    This relies on a rational approximation (error ~ 2e-7 (x<1) or 5e-5 (x>1)
207#
208proc ::math::special::exponential_E1 {x} {
209    if { $x <= 0.0 } {
210        error "Domain error: x must be positive"
211    }
212
213   if { $x < 1.0 } {
214      return [expr {-log($x)+((((0.00107857*$x-0.00976004)*$x+0.05519968)*$x-0.24991055)*$x+0.99999193)*$x-0.57721566}]
215   } else {
216      set xexpe [expr {($x*$x+2.334733*$x+0.250621)/($x*$x+3.330657*$x+1.681534)}]
217      return [expr {$xexpe/($x*exp($x))}]
218   }
219}
220
221# exponential_li --
222#    Compute the logarithmic integral
223# Arguments:
224#    x       Value of the argument
225# Result:
226#    Value of the integral 1/ln(x) from 0 to x
227#
228proc ::math::special::exponential_li {x} {
229    if { $x < 0 } {
230        return -code error "Argument must be positive or zero"
231    } else {
232        if { $x == 0.0 } {
233            return 0.0
234        } else {
235            return [exponential_Ei [expr {log($x)}]]
236        }
237    }
238}
239
240# exponential_Shi --
241#    Compute the hyperbolic sine integral
242# Arguments:
243#    x       Value of the argument
244# Result:
245#    Value of the integral sinh(x)/x from 0 to x
246#
247proc ::math::special::exponential_Shi {x} {
248    if { $x < 0 } {
249        return -code error "Argument must be positive or zero"
250    } else {
251        if { $x == 0.0 } {
252            return 0.0
253        } else {
254            proc g {x} {
255               return [expr {sinh($x)/$x}]
256            }
257            return [lindex [::math::calculus::romberg g 0.0 $x] 0]
258        }
259    }
260}
261
262# exponential_Chi --
263#    Compute the hyperbolic cosine integral
264# Arguments:
265#    x       Value of the argument
266# Result:
267#    Value of the integral (cosh(x)-1)/x from 0 to x
268#
269proc ::math::special::exponential_Chi {x} {
270    variable gamma
271    if { $x < 0 } {
272        return -code error "Argument must be positive or zero"
273    } else {
274        if { $x == 0.0 } {
275            return 0.0
276        } else {
277            proc g {x} {
278               return [expr {(cosh($x)-1.0)/$x}]
279            }
280            set integral [lindex [::math::calculus::romberg g 0.0 $x] 0]
281            return [expr {$gamma+log($x)+$integral}]
282        }
283    }
284}
285
286# exponential_Si --
287#    Compute the sine integral
288# Arguments:
289#    x       Value of the argument
290# Result:
291#    Value of the integral sin(x)/x from 0 to x
292#
293proc ::math::special::exponential_Si {x} {
294    variable halfpi
295    if { $x < 0 } {
296        return -code error "Argument must be positive or zero"
297    } else {
298        if { $x == 0.0 } {
299            return 0.0
300        } else {
301            if { $x < 1.0 } {
302                proc g {x} {
303                    return [expr {sin($x)/$x}]
304                }
305                return [lindex [::math::calculus::romberg g 0.0 $x] 0]
306            } else {
307                foreach {f g} [ComputeExponFG $x] {break}
308                return [expr {$halfpi-$f*cos($x)-$g*sin($x)}]
309            }
310        }
311    }
312}
313
314# exponential_Ci --
315#    Compute the cosine integral
316# Arguments:
317#    x       Value of the argument
318# Result:
319#    Value of the integral (cosh(x)-1)/x from 0 to x
320#
321proc ::math::special::exponential_Ci {x} {
322    variable gamma
323
324    if { $x < 0 } {
325        return -code error "Argument must be positive or zero"
326    } else {
327        if { $x == 0.0 } {
328            return 0.0
329        } else {
330            if { $x < 1.0 } {
331                proc g {x} {
332                    return [expr {(cos($x)-1.0)/$x}]
333                }
334                set integral [lindex [::math::calculus::romberg g 0.0 $x] 0]
335                return [expr {$gamma+log($x)+$integral}]
336            } else {
337                foreach {f g} [ComputeExponFG $x] {break}
338                return [expr {$f*sin($x)-$g*cos($x)}]
339            }
340        }
341    }
342}
343
344# some tests --
345#    Reproduce tables 5.1, 5.2 from Abramowitz & Stegun,
346
347if { [info exists ::argv0] && ![string compare $::argv0 [info script]] } {
348namespace eval ::math::special {
349for { set i 0.01 } { $i < 0.505 } { set i [expr { $i + 0.01 }] } {
350    set ei [exponential_Ei $i]
351    set e1 [expr { - [exponential_Ei [expr { - $i }]] }]
352    puts [format "%9.6f\t%.10g\t%.10g" $i \
353              [expr {($ei - log($i) - 0.57721566490153286)/$i} ] \
354              [expr {($e1 + log($i) + 0.57721566490153286) / $i }]]
355}
356puts {}
357for { set i 0.5 } { $i < 2.005 } { set i [expr { $i + 0.01 }] } {
358    set ei [exponential_Ei $i]
359    set e1 [expr { - [exponential_Ei [expr { - $i }]] }]
360    puts [format "%9.6f\t%.10g\t%.10g" $i $ei $e1]
361}
362puts {}
363for {} { $i < 10.05 } { set i [expr { $i + 0.1 }] } {
364    set ei [exponential_Ei $i]
365    set e1 [expr { - [exponential_Ei [expr { - $i }]] }]
366    puts [format "%9.6f\t%.10g\t%.10g" $i \
367              [expr { $i * exp(-$i) * $ei }] \
368              [expr { $i * exp($i) * $e1 }]]
369
370}
371puts {}
372for {set ooi 0.1} { $ooi > 0.0046 } { set ooi [expr { $ooi - 0.005 }] } {
373    set i [expr { 1.0 / $ooi }]
374    set ri [expr { round($i) }]
375    set ei [exponential_Ei $i]
376    set e1 [expr { - [exponential_Ei [expr { - $i }]] }]
377    puts [format "%9.6f\t%.10g\t%.10g\t%d" $i \
378              [expr { $i * exp(-$i) * $ei }] \
379              [expr { $i * exp($i) * $e1 }] \
380              $ri]
381}
382puts {}
383
384# Reproduce table 5.4 from Abramowitz and Stegun
385
386for { set x 0.00 } { $x < 0.505 } { set x [expr { $x + 0.01 }] } {
387    set line [format %4.2f $x]
388    if { $x == 0.00 } {
389        append line { } 1.0000000
390    } else {
391        append line { } [format %9.7f \
392                             [expr { [exponential_En 2 $x] - $x * log($x) }]]
393    }
394    foreach n { 3 4 10 20 } {
395        append line { } [format %9.7f [exponential_En $n $x]]
396    }
397    puts $line
398}
399puts {}
400for { set x 0.50 } { $x < 2.005 } { set x [expr { $x + 0.01 }] } {
401    set line [format %4.2f $x]
402    foreach n { 2 3 4 10 20 } {
403        append line { } [format %9.7f [exponential_En $n $x]]
404    }
405    puts $line
406}
407puts {}
408
409for { set oox 0.5 } { $oox > 0.1025 } { set oox [expr { $oox - 0.05 }] } {
410    set line [format %4.2f $oox]
411    set x [expr { 1.0 / $oox }]
412    set rx [expr { round( $x ) }]
413    foreach n { 2 3 4 10 20 } {
414        set en [exponential_En $n [expr { 1.0 / $oox }]]
415        append line { } [format %9.7f [expr { ( $x + $n ) * exp($x) * $en }]]
416    }
417    append line { } [format %3d $rx]
418    puts $line
419}
420for { set oox 0.10 } { $oox > 0.005 } { set oox [expr { $oox - 0.01 }] } {
421    set line [format %4.2f $oox]
422    set x [expr { 1.0 / $oox }]
423    set rx [expr { round( $x ) }]
424    foreach n { 2 3 4 10 20 } {
425        set en [exponential_En $n $x]
426        append line { } [format %9.7f [expr { ( $x + $n ) * exp($x) * $en }]]
427    }
428    append line { } [format %3d $rx]
429    puts $line
430}
431puts {}
432catch {exponential_Ei 0.0} result; puts $result
433}
434}
435