1# special.tcl -- 2# Provide well-known special mathematical functions 3# 4# This file contains a collection of tests for one or more of the Tcllib 5# procedures. Sourcing this file into Tcl runs the tests and 6# generates output for errors. No output means no errors were found. 7# 8# Copyright (c) 2004 by Arjen Markus. All rights reserved. 9# 10# RCS: @(#) $Id: special.tcl,v 1.13 2008/08/13 07:28:47 arjenmarkus Exp $ 11# 12package require math 13package require math::constants 14package require math::statistics 15 16# namespace special 17# Create a convenient namespace for the "special" mathematical functions 18# 19namespace eval ::math::special { 20 # 21 # Define a number of common mathematical constants 22 # 23 ::math::constants::constants pi 24 variable halfpi [expr {$pi/2.0}] 25 26 # 27 # Functions defined in other math submodules 28 # 29 if { [info commands Beta] == {} } { 30 namespace import ::math::Beta 31 namespace import ::math::ln_Gamma 32 } 33 34 # 35 # Export the various functions 36 # 37 namespace export Beta ln_Gamma Gamma erf erfc fresnel_C fresnel_S sinc 38} 39 40# Gamma -- 41# The Gamma function - synonym for "factorial" 42# 43proc ::math::special::Gamma {x} { 44 if { [catch { expr {exp( [ln_Gamma $x] )} } result] } { 45 return -code error -errorcode $::errorCode $result 46 } 47 return $result 48} 49 50# erf -- 51# The error function 52# Arguments: 53# x The value for which the function must be evaluated 54# Result: 55# erf(x) 56# Note: 57# The algoritm used is due to George Marsaglia 58# See: http://www.velocityreviews.com/forums/t317358-erf-function-in-c.html 59# I did not want to copy and convert the even more accurate but 60# rather lengthy algorithm used by lcc-win32/Sun 61# 62proc ::math::special::erf {x} { 63 set x [expr {$x*sqrt(2.0)}] 64 65 if { $x > 10.0 } { return 1.0 } 66 if { $x < -10.0 } { return -1.0 } 67 68 set a 1.2533141373155 69 set b -1.0 70 set pwr 1.0 71 set t 0.0 72 set z 0.0 73 74 set s [expr {$a+$b*$x}] 75 76 set i 2 77 while { $s != $t } { 78 set a [expr {($a+$z*$b)/double($i)}] 79 set b [expr {($b+$z*$a)/double($i+1)}] 80 set pwr [expr {$pwr*$x*$x}] 81 set t $s 82 set s [expr {$s+$pwr*($a+$x*$b)}] 83 84 incr i 2 85 } 86 87 return [expr {1.0-2.0*$s*exp(-0.5*$x*$x-0.9189385332046727418)}] 88} 89 90 91 92# erfc -- 93# The complement of the error function 94# Arguments: 95# x The value for which the function must be evaluated 96# Result: 97# erfc(x) = 1.0-erf(x) 98# 99proc ::math::special::erfc {x} { 100 set x [expr {$x*sqrt(2.0)}] 101 102 if { $x > 10.0 } { return 0.0 } 103 if { $x < -10.0 } { return 0.0 } 104 105 set a 1.2533141373155 106 set b -1.0 107 set pwr 1.0 108 set t 0.0 109 set z 0.0 110 111 set s [expr {$a+$b*$x}] 112 113 set i 2 114 while { $s != $t } { 115 set a [expr {($a+$z*$b)/double($i)}] 116 set b [expr {($b+$z*$a)/double($i+1)}] 117 set pwr [expr {$pwr*$x*$x}] 118 set t $s 119 set s [expr {$s+$pwr*($a+$x*$b)}] 120 121 incr i 2 122 } 123 124 return [expr {2.0*$s*exp(-0.5*$x*$x-0.9189385332046727418)}] 125} 126 127 128# ComputeFG -- 129# Compute the auxiliary functions f and g 130# 131# Arguments: 132# x Parameter of the integral (x>=0) 133# Result: 134# Approximate values for f and g 135# Note: 136# See Abramowitz and Stegun. The accuracy is 2.0e-3. 137# 138proc ::math::special::ComputeFG {x} { 139 list [expr {(1.0+0.926*$x)/(2.0+1.792*$x+3.104*$x*$x)}] \ 140 [expr {1.0/(2.0+4.142*$x+3.492*$x*$x+6.670*$x*$x*$x)}] 141} 142 143# fresnel_C -- 144# Compute the Fresnel cosine integral 145# 146# Arguments: 147# x Parameter of the integral (x>=0) 148# Result: 149# Value of C(x) = integral from 0 to x of cos(0.5*pi*x^2) 150# Note: 151# This relies on a rational approximation of the two auxiliary functions f and g 152# 153proc ::math::special::fresnel_C {x} { 154 variable halfpi 155 if { $x < 0.0 } { 156 error "Domain error: x must be non-negative" 157 } 158 159 if { $x == 0.0 } { 160 return 0.0 161 } 162 163 foreach {f g} [ComputeFG $x] {break} 164 165 set xarg [expr {$halfpi*$x*$x}] 166 167 return [expr {0.5+$f*sin($xarg)-$g*cos($xarg)}] 168} 169 170# fresnel_S -- 171# Compute the Fresnel sine integral 172# 173# Arguments: 174# x Parameter of the integral (x>=0) 175# Result: 176# Value of S(x) = integral from 0 to x of sin(0.5*pi*x^2) 177# Note: 178# This relies on a rational approximation of the two auxiliary functions f and g 179# 180proc ::math::special::fresnel_S {x} { 181 variable halfpi 182 if { $x < 0.0 } { 183 error "Domain error: x must be non-negative" 184 } 185 186 if { $x == 0.0 } { 187 return 0.0 188 } 189 190 foreach {f g} [ComputeFG $x] {break} 191 192 set xarg [expr {$halfpi*$x*$x}] 193 194 return [expr {0.5-$f*cos($xarg)-$g*sin($xarg)}] 195} 196 197# sinc -- 198# Compute the sinc function 199# Arguments: 200# x Value of the argument 201# Result: 202# sin(x)/x 203# 204proc ::math::special::sinc {x} { 205 if { $x == 0.0 } { 206 return 1.0 207 } else { 208 return [expr {sin($x)/$x}] 209 } 210} 211 212# Bessel functions and elliptic integrals -- 213# 214source [file join [file dirname [info script]] "bessel.tcl"] 215source [file join [file dirname [info script]] "classic_polyns.tcl"] 216source [file join [file dirname [info script]] "elliptic.tcl"] 217source [file join [file dirname [info script]] "exponential.tcl"] 218 219package provide math::special 0.2.2 220