[Arjen Markus] (22 november 2002) It must be the time of year, but I felt very mathematical this week. So I sat down and wrote a little script that does a mathematical job: it allows you to define polynomial functions and (truncated) power series. Obviously, this is just a quick version and a full-blown package should include such things as:

   * Adding, multiplying polynomials (division with remainder)
   * Differentiation, integration (the result being another polynomial function)
   * Scaling
   * Introspection and error estimates 
   * Numerous other fascinating operations

Nevertheless, if you have a need for other functions than the ones offered by [expr], this may give you some idea. The proof of the pudding is usually in the eating - so the example includes the zeroth order Bessel function of the first kind, J0. It is my favourite :)

[AM] A more direct evaluation of Bessel functions can be found at: [Bessel functions]

----
 # polynoms.tcl -- 
 # 
 #    Package for handling polynoms and power series
 #    (sample Workbench module)
 #
 # Notes:
 #    This package is a quick hack to get started only
 #
 # Version information:
 #    version 0.1: initial implementation, november 2002

 package provide Polynoms 0.1

 namespace eval ::polynoms {

    namespace export polynom powerseries

 # powerseries --
 #    Construct a procedure that implements a power series
 #    and return its name
 #
 # Arguments:
 #    name         Name of the new procedure
 #    nterms       Number of terms
 #    calcterm     Procedure to calculate a term
 #
 # Result:
 #    Name of procedure that will calculate the power series
 #
 proc powerseries {name nterms calcterm} {

    set terms {}
    for { set n 0 } { $n <= $nterms } { incr n } {
       lappend terms [$calcterm $n]
    }

    polynom $name $terms
 }

 # polynom --
 #    Construct a procedure that implements a polynom
 #    and return its name
 #
 # Arguments:
 #    name         Name of the new procedure
 #    terms        Coefficients (from degree 0 to degree n)
 #
 # Result:
 #    Name of procedure that will calculate the polygon
 #
 proc polynom {name terms} {

    set first  [lindex $terms end]
    set rterms {}
    for { set n [expr {[llength $terms]-2}] } { $n >= 0 } { incr n -1 } {
       lappend rterms [lindex $terms $n]
    }

    interp alias {} $name {} [namespace current]::EvalPolynom \
       $first $rterms
    return $name
 }

 # EvalPolynom --
 #    Calculate the value of a polynom
 #
 # Arguments:
 #    rterms       Coefficients from degree n to degree 0
 #    x            Value of the argument
 #
 # Result:
 #    Value of the polynom at x
 #
 proc EvalPolynom {first rterms x} {
    set value $first
    foreach term $rterms {
       set value [expr {$value*$x+$term}]
    }
    return $value
 }
 
 } ;# End of namespace
 
 
 #
 # Run the program
 #
 namespace import ::polynoms::*

 puts "Simple parabola"
 polynom parabola {0.0 0.0 1.0}
 for { set i 0 } { $i < 10 } { incr i } {
    puts [parabola $i]
 }

 puts "Exponential function"
 
 proc exp_terms {n} {
    set fac 1
    for {set i 1} {$i <= $n} {incr i} {
       set fac [expr {$fac*$i}]
    }
    expr {1.0/$fac}
 } 
 powerseries exp 10 exp_terms
 for { set i 0 } { $i < 10 } { incr i } {
    puts "[exp $i] [expr {exp($i)}]"
 }
 
 puts "Bessel function J0"

 proc J0_terms {n} {
    if { $n%2 == 1 } {
       return 0.0
    } else {
       set fac 1
       set sign [expr {$n%4==2?-1:1}]
       for {set i 1} {$i <= $n/2} {incr i} {
          set fac  [expr {$fac*$i}]
          set sign [expr {$sign/4.0}]
       }
       expr {$sign/($fac*$fac)}
    }
 }
 powerseries J0 10 J0_terms
 for { set x 0.0 } { $x < 2.0 } { set x [expr {$x+0.2}] } {
    puts "$x [J0 $x]"
 }

----
[AM] (13 july 2004) A more mature package is on its way for Tcllib. See [Evaluating polynomial functions] for some numerical experiments.
----
[Category Mathematics]