1# classic_polyns.tcl --
2#    Implement procedures for the classic orthogonal polynomials
3#
4package require math::polynomials
5
6namespace eval ::math::special {
7    if {[info commands addPolyn] == {} } {
8        namespace import ::math::polynomials::*
9    }
10}
11
12
13# legendre --
14#    Return the nth degree Legendre polynomial
15#
16# Arguments:
17#    n            The degree of the polynomial
18# Result:
19#    Polynomial definition
20#
21proc ::math::special::legendre {n} {
22    if { ! [string is integer -strict $n] || $n < 0 } {
23        return -code error "Degree must be a non-negative integer"
24    }
25
26    set pnm1   [polynomial 1.0]
27    set pn     [polynomial {0.0 1.0}]
28
29    if { $n == 0 } {return $pnm1}
30    if { $n == 1 } {return $pn}
31
32    set degree 1
33    while { $degree < $n } {
34        set an       [expr {(2.0*$degree+1.0)/($degree+1.0)}]
35        set bn       0.0
36        set cn       [expr {$degree/($degree+1.0)}]
37        set factor_n [polynomial [list $bn $an]]
38        set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]]
39        set term_n   [multPolyn $factor_n $pn]
40        set pnp1     [addPolyn $term_n $term_nm1]
41
42        set pnm1     $pn
43        set pn       $pnp1
44        incr degree
45    }
46
47    return $pnp1
48}
49
50# chebyshev --
51#    Return the nth degree Chebeyshev polynomial of the first kind
52#
53# Arguments:
54#    n            The degree of the polynomial
55# Result:
56#    Polynomial definition
57#
58proc ::math::special::chebyshev {n} {
59    if { ! [string is integer -strict $n] || $n < 0 } {
60        return -code error "Degree must be a non-negative integer"
61    }
62
63    set pnm1   [polynomial 1.0]
64    set pn     [polynomial {0.0 1.0}]
65
66    if { $n == 0 } {return $pnm1}
67    if { $n == 1 } {return $pn}
68
69    set degree 1
70    while { $degree < $n } {
71        set an       2.0
72        set bn       0.0
73        set cn       1.0
74        set factor_n [polynomial [list $bn $an]]
75        set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]]
76        set term_n   [multPolyn $factor_n $pn]
77        set pnp1     [addPolyn $term_n $term_nm1]
78
79        set pnm1     $pn
80        set pn       $pnp1
81        incr degree
82    }
83
84    return $pnp1
85}
86
87# laguerre --
88#    Return the nth degree Laguerre polynomial with parameter alpha
89#
90# Arguments:
91#    alpha        The parameter for the polynomial
92#    n            The degree of the polynomial
93# Result:
94#    Polynomial definition
95#
96proc ::math::special::laguerre {alpha n} {
97    if { ! [string is double -strict $alpha] } {
98        return -code error "Parameter must be a double"
99    }
100    if { ! [string is integer -strict $n] || $n < 0 } {
101        return -code error "Degree must be a non-negative integer"
102    }
103
104    set pnm1   [polynomial 1.0]
105    set pn     [polynomial [list [expr {1.0-$alpha}] -1.0]]
106
107    if { $n == 0 } {return $pnm1}
108    if { $n == 1 } {return $pn}
109
110    set degree 1
111    while { $degree < $n } {
112        set an       [expr {-1.0/($degree+1.0)}]
113        set bn       [expr {(2.0*$degree+$alpha+1)/($degree+1.0)}]
114        set cn       [expr {($degree+$alpha)/($degree+1.0)}]
115        set factor_n [polynomial [list $bn $an]]
116        set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]]
117        set term_n   [multPolyn $factor_n $pn]
118        set pnp1     [addPolyn $term_n $term_nm1]
119
120        set pnm1     $pn
121        set pn       $pnp1
122        incr degree
123    }
124
125    return $pnp1
126}
127
128# hermite --
129#    Return the nth degree Hermite polynomial
130#
131# Arguments:
132#    n            The degree of the polynomial
133# Result:
134#    Polynomial definition
135#
136proc ::math::special::hermite {n} {
137    if { ! [string is integer -strict $n] || $n < 0 } {
138        return -code error "Degree must be a non-negative integer"
139    }
140
141    set pnm1   [polynomial 1.0]
142    set pn     [polynomial {0.0 2.0}]
143
144    if { $n == 0 } {return $pnm1}
145    if { $n == 1 } {return $pn}
146
147    set degree 1
148    while { $degree < $n } {
149        set an       2.0
150        set bn       0.0
151        set cn       [expr {2.0*$degree}]
152        set factor_n [polynomial [list $bn $an]]
153        set term_n   [multPolyn $factor_n $pn]
154        set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]]
155        set pnp1     [addPolyn $term_n $term_nm1]
156
157        set pnm1     $pn
158        set pn       $pnp1
159        incr degree
160    }
161
162    return $pnp1
163}
164
165# some tests --
166#
167if { 0 } {
168set tcl_precision 17
169
170puts "Legendre:"
171foreach n {0 1 2 3 4} {
172    puts [::math::special::legendre $n]
173}
174
175puts "Chebyshev:"
176foreach n {0 1 2 3 4} {
177    puts [::math::special::chebyshev $n]
178}
179
180puts "Laguerre (alpha=0):"
181foreach n {0 1 2 3 4} {
182    puts [::math::special::laguerre 0.0 $n]
183}
184puts "Laguerre (alpha=1):"
185foreach n {0 1 2 3 4} {
186    puts [::math::special::laguerre 1.0 $n]
187}
188
189puts "Hermite:"
190foreach n {0 1 2 3 4} {
191    puts [::math::special::hermite $n]
192}
193}
194