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