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