Recursively calculate Chebyshev polynomials

chebyshev.tcl

#!/usr/bin/tclsh

# we need Tcl 8.5 for lrepeat to be defined
package require Tcl 8.5

#
# recursively calculate Chebyshev polynomials
#
#   (see http://en.wikipedia.org/wiki/Chebyshev_polynomials)
#

proc min {x y} {
    if {$x <= $y} {return $x} {return $y}
}

#
# NOTE: polynomials are represented as lists, e.g. the polynomial 2x^3 + 3x + 1
#       is represented as [1, 3, 0, 2]
#

#
# trim out zero term of higher degrees, e.g. 0x^4 + x^3 + 1 => x^3 + 1
#
proc poly_trim {p} {
    for {set i [expr {[llength $p] - 1}]} {$i >= 0} {set i [expr {$i - 1}]} {
        if {[lindex $p $i] == 0} {
            set p [lreplace $p end end]
        } else {
            break
        }
    }

    return $p
}

#
# add two polynomials
#
proc poly_add {p q} {
    set m [min [llength $p] [llength $q]]
    set sum [list]

    for {set i 0} {$i < $m} {incr i} {
        lappend sum [expr {[lindex $p $i] + [lindex $q $i]}]
    }

    if {[llength $p] > [llength $q]} {
        for {set i $m} {$i < [llength $p]} {incr i} {
            lappend sum [lindex $p $i]
        }
    } elseif {[llength $p] < [llength $q]} {
        for {set i $m} {$i < [llength $q]} {incr i} {
            lappend sum [lindex $q $i]
        }
    }

    return [poly_trim $sum]
}

#
# multiply two polynomials
#
proc poly_multiply {p q} {
    set product [lrepeat [expr {[llength $p] + [llength $q] - 1}] 0]

    for {set i 0} {$i < [llength $p]} {incr i} {
        for {set j 0} {$j < [llength $q]} {incr j} {
            set k [expr {$i + $j}]

            set product [
                lreplace $product $k $k [
                    expr {
                        [lindex $product $k]
                        +
                        [expr {[lindex $p $i] * [lindex $q $j]}]
                    }
                ]
            ]
        }
    }

    return [poly_trim $product]
}

#
# Chebyshev polynomials of the first kind
#
# these are defined recursively as follows:
#
#   T_0(x) = 1
#   T_1(x) = x
#   T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x)
#
proc chebyshev1 {n} {
    if {$n <= 0} {
        return [list 1]
    } elseif {$n == 1} {
        return [list 0 1]
    } else {
        return [
            poly_add [
                poly_multiply [list 0 2] [chebyshev1 [expr {$n - 1}]]
            ] [
                poly_multiply [list -1] [chebyshev1 [expr {$n - 2}]]
            ]
        ]
    }
}

#
# Chebyshev polynomials of the second kind
#
# these are defined recursively as follows:
#
#   U_0(x) = 1
#   U_1(x) = 2x
#   U_{n+1}(x) = 2xU_n(x) - U_{n-1}(x)
#
proc chebyshev2 {n} {
    if {$n <= 0} {
        return [list 1]
    } elseif {$n == 1} {
        return [list 0 2]
    } else {
        return [
            poly_add [
                poly_multiply [list 0 2] [chebyshev2 [expr {$n - 1}]]
            ] [
                poly_multiply [list -1] [chebyshev2 [expr {$n - 2}]]
            ]
        ]
    }
}


puts [chebyshev1 10]
puts [chebyshev2 10]

# ex: ts=4 sw=4 et filetype=tcl noexpandtab