1 files changed, 200 insertions, 0 deletions
diff --git a/tcllib/modules/math/classic_polyns.tcl b/tcllib/modules/math/classic_polyns.tcl
new file mode 100755
index 0000000..1fd9cd0
--- /dev/null
+++ b/tcllib/modules/math/classic_polyns.tcl
@@ -0,0 +1,200 @@
+# classic_polyns.tcl --
+#    Implement procedures for the classic orthogonal polynomials
+#
+package require math::polynomials
+
+namespace eval ::math::special {
+    if {[info commands addPolyn] == {} } {
+        namespace import ::math::polynomials::*
+    }
+}
+
+
+# legendre --
+#    Return the nth degree Legendre polynomial
+#
+# Arguments:
+#    n            The degree of the polynomial
+# Result:
+#    Polynomial definition
+#
+proc ::math::special::legendre {n} {
+    if { ! [string is integer -strict $n] || $n < 0 } {
+        return -code error "Degree must be a non-negative integer"
+    }
+
+    set pnm1   [polynomial 1.0]
+    set pn     [polynomial {0.0 1.0}]
+
+    if { $n == 0 } {return $pnm1}
+    if { $n == 1 } {return $pn}
+
+    set degree 1
+    while { $degree < $n } {
+        set an       [expr {(2.0*$degree+1.0)/($degree+1.0)}]
+        set bn       0.0
+        set cn       [expr {$degree/($degree+1.0)}]
+        set factor_n [polynomial [list $bn $an]]
+        set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]]
+        set term_n   [multPolyn $factor_n $pn]
+        set pnp1     [addPolyn $term_n $term_nm1]
+
+        set pnm1     $pn
+        set pn       $pnp1
+        incr degree
+    }
+
+    return $pnp1
+}
+
+# chebyshev --
+#    Return the nth degree Chebeyshev polynomial of the first kind
+#
+# Arguments:
+#    n            The degree of the polynomial
+# Result:
+#    Polynomial definition
+#
+proc ::math::special::chebyshev {n} {
+    if { ! [string is integer -strict $n] || $n < 0 } {
+        return -code error "Degree must be a non-negative integer"
+    }
+
+    set pnm1   [polynomial 1.0]
+    set pn     [polynomial {0.0 1.0}]
+
+    if { $n == 0 } {return $pnm1}
+    if { $n == 1 } {return $pn}
+
+    set degree 1
+    while { $degree < $n } {
+        set an       2.0
+        set bn       0.0
+        set cn       1.0
+        set factor_n [polynomial [list $bn $an]]
+        set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]]
+        set term_n   [multPolyn $factor_n $pn]
+        set pnp1     [addPolyn $term_n $term_nm1]
+
+        set pnm1     $pn
+        set pn       $pnp1
+        incr degree
+    }
+
+    return $pnp1
+}
+
+# laguerre --
+#    Return the nth degree Laguerre polynomial with parameter alpha
+#
+# Arguments:
+#    alpha        The parameter for the polynomial
+#    n            The degree of the polynomial
+# Result:
+#    Polynomial definition
+#
+proc ::math::special::laguerre {alpha n} {
+    if { ! [string is double -strict $alpha] } {
+        return -code error "Parameter must be a double"
+    }
+    if { ! [string is integer -strict $n] || $n < 0 } {
+        return -code error "Degree must be a non-negative integer"
+    }
+
+    set pnm1   [polynomial 1.0]
+    set pn     [polynomial [list [expr {1.0-$alpha}] -1.0]]
+
+    if { $n == 0 } {return $pnm1}
+    if { $n == 1 } {return $pn}
+
+    set degree 1
+    while { $degree < $n } {
+        set an       [expr {-1.0/($degree+1.0)}]
+        set bn       [expr {(2.0*$degree+$alpha+1)/($degree+1.0)}]
+        set cn       [expr {($degree+$alpha)/($degree+1.0)}]
+        set factor_n [polynomial [list $bn $an]]
+        set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]]
+        set term_n   [multPolyn $factor_n $pn]
+        set pnp1     [addPolyn $term_n $term_nm1]
+
+        set pnm1     $pn
+        set pn       $pnp1
+        incr degree
+    }
+
+    return $pnp1
+}
+
+# hermite --
+#    Return the nth degree Hermite polynomial
+#
+# Arguments:
+#    n            The degree of the polynomial
+# Result:
+#    Polynomial definition
+#
+proc ::math::special::hermite {n} {
+    if { ! [string is integer -strict $n] || $n < 0 } {
+        return -code error "Degree must be a non-negative integer"
+    }
+
+    set pnm1   [polynomial 1.0]
+    set pn     [polynomial {0.0 2.0}]
+
+    if { $n == 0 } {return $pnm1}
+    if { $n == 1 } {return $pn}
+
+    set degree 1
+    while { $degree < $n } {
+        set an       2.0
+        set bn       0.0
+        set cn       [expr {2.0*$degree}]
+        set factor_n [polynomial [list $bn $an]]
+        set term_n   [multPolyn $factor_n $pn]
+        set term_nm1 [multPolyn $pnm1 [expr {-1.0*$cn}]]
+        set pnp1     [addPolyn $term_n $term_nm1]
+
+        set pnm1     $pn
+        set pn       $pnp1
+        incr degree
+    }
+
+    return $pnp1
+}
+
+# some tests --
+#
+if { 0 } {
+set prec $::tcl_precision
+if {![package vsatisfies [package provide Tcl] 8.5]} {
+    set ::tcl_precision 17
+} else {
+    set ::tcl_precision 0
+}
+
+puts "Legendre:"
+foreach n {0 1 2 3 4} {
+    puts [::math::special::legendre $n]
+}
+
+puts "Chebyshev:"
+foreach n {0 1 2 3 4} {
+    puts [::math::special::chebyshev $n]
+}
+
+puts "Laguerre (alpha=0):"
+foreach n {0 1 2 3 4} {
+    puts [::math::special::laguerre 0.0 $n]
+}
+puts "Laguerre (alpha=1):"
+foreach n {0 1 2 3 4} {
+    puts [::math::special::laguerre 1.0 $n]
+}
+
+puts "Hermite:"
+foreach n {0 1 2 3 4} {
+    puts [::math::special::hermite $n]
+}
+
+set ::tcl_precision $prec
+}