Mathematics with Tcl

Arjen Markus (20 december 2007) I wrote this page for the website Torsten Berg is designing

Lars H, 2007-12-21: Do you accept changes, or would suggestions be better? Some points that caught my eye:

A group action is when you let a group act on some set. For composition of permutations I would rather speak of group operation (what you find out is roughly that the two are the same for this representation of permutations).
It could be mentioned that the taxicab solution includes the most common approach to sets in Tcl (and one that can even be extended to multisets).
Shouldn't the last line of the diff example results be 200 40000 - 399.999997171?
The examples chosen very often end up in automatically generating code, which (IMHO) is not that typical for mathematics in Tcl, and thus could give a wrong first impression (e.g. "the core language is not suited for mathematics, so you have to do a bunch of metaprogramming instead"). The implementation in memoizing would probably show that idea better.

AM (21 december 2007) Sure, I accept changes and suggestions ;). I had not noticed that I generate a lot of automatic code. Perhaps it is time for a review. I will certainly have a closer look! And changes are welcome.

Mathematics with Tcl/Tk

This is an introduction to Tcl/Tk with an emphasis on mathematics. It is assumed that you are more or less familiar with the basics of Tcl and at the very least with programming in general - variables, loops, functions, if-statements, they should all be familiar for you.

We will look at both numerical analysis and at such more abstract topics like group theory. Rather than develop a fully functional package for any of the topics, we will simply show how you can use Tcl and the standard data structures it has to solve mathematical problems.

Mathematics is simply such a vast subject, it is impossible to do more than look at a few things within the space of a tutorial. The emphasis will be on data structures that can help solve the particular problem we will look at.

Note: We will assume Tcl 8.4 in this introduction with a number of packages. Some remarks will be made about Tcl 8.5, as that offers even more facilities.

Group theory: permutations

Group theory is a very important part of modern algebra. And within group theory permutations take a special place, as the group of all permutations of N objects has as its subgroups all possible groups of order N. So permutations seem like a nice starting point for looking at group theory.

The first thing we need is a way to represent a permutation. This is very naturally done with a list that contains the new positions for all objects:

Suppose permutation P of five objects sends the first object to the second position, the second object to the fifth, third remains in the same position, the fourth is sent to the first position and finally the fifth is put on the fourth position:

            P
   ABCDE ------> DACEB

We can use the list {2 5 3 1 4} to represent this permutation. Because Tcl uses 0 as the first index for list items a more convenient way to represent P is: {1 4 2 0 3}

The following procedure will reorder a list of five things according to the given permutation:

    proc perm {permute objects} {
        set result $objects

        foreach p $permute o $objects {
            lset result $p $o
        }
        return $result
    }

So:

    % puts [perm {1 4 2 0 3} {A B C D E}
    ==> D A C E B

The next step is to see how we can combine permutations - that is: we want to implement the group action.

Actually, that is trivial: to get a permutation that is the combined effect of permutation A and B (first permute by B, then by A), just use the procedure perm:

    set AB [perm $A $B]

The law of associativity is expressed as follows:

    set ABC1 [perm [perm $A $B] $C]
    set ABC2 [perm $A [perm $B $C]]

    if { $ABC1 != $ABC2 } {
        puts "Violation of the law of associativity!"
    }

Just for fun, let us compute the subgroup that arises with the permutations {0 1 4 2 3} and {2 1 0 3 4}:

    proc computeGroup {operation elements} {
        set newelems 0

        #
        # We want unique elements. Use an array to do so:
        # the name of the array element is automatically unique
        #
        foreach e $elements {
            set result($e) 1
        }

        #
        # We multiply each element with each other element
        # and store the result
        #
        foreach a $elements {
            foreach b $elements {
                set c [$operation $a $b]
                if { ! [info exists result($c)] } {
                    set newelems   1
                    set result($c) 1
                }
            }
        }

        #
        # Now we have a new array, possibly with new elements
        # If so, restart the procedure. If not, return a list of
        # elements
        #
        set elements [array names result]

        if { $newelems } {
            return [computeGroup $operation $elements]
        } else {
            return $elements
        }
    }

So compute the subgroup we are looking for:

    set group [computeGroup perm {{1 0 4 2 3} {2 1 0 3 4}}]

and the result:

(yes, it is a group of 24 elements, actually the group of all permutations of four elements, as the element in the second position is not moved and all others are).

We could have implemented the computeGroup in other ways as well, for instance:

Use a list to collect all elements and only add an element if it does not occur yet
Use a list to collect all elements, sort the list retaining unique elements only

but the property of Tcl arrays to have only one element of a specific name makes them rather suitable for this kind of book keeping.

Note: Tcl 8.5 introduces dictionaries, these are even more suitable for this type of work, as they can be passed as arguments, whereas arrays can only be passed by name.

NOTE: use dictionaries to implement sets?

Taxicab numbers

(See [L1 ] for an explanation of the name - something I didn't know until coincidentally reading about it a couple of days ago - CJL)

Taxicab numbers are numbers that can be represented in at least two ways as the sum of positive cubes. There is another definition too, but let us focus on this one. The smallest such number is 1729:

    1729 = 1 + 12^3 = 10^3 + 9^3

Finding them is just a matter of brute force:

Compute x^3+y^3 where x goes from 1 to some limit and y from 1 to x (no need to go further)
Then scan the list of resulting sums for values that appear more than once

We could do it literally in this way ... But what about the following method:

We use an associative array (Tcl's arrays are associative) with the sum x^3+y^3 as the key
We add the pairs x,y to the value of each key
When we have worked our way through the double loop, we check for array elements that have two or more pairs
No need to check the whole list for any doubles, that is taken care of by the key automatically

Here is the code:

   proc taxicab {limit} {
       for { set x 1 } { $x <= $limit } { incr x } {
           for { set y 1 } { $y <= $x } { incr y } {
               set z [expr {$x*$x*$x+$y*$y*$y}]
               lappend taxicab($z) [list $x $y]
           }
       }

       foreach {z xypairs} [array get taxicab] {
           if { [llength $xypairs] > 1 } {
               puts "$z: [join $xypairs " | "]"
           }
       }
   }

   taxicab 100

The output this produces is:

    216125: 50 45 | 60 5
    558441: 72 57 | 81 30
    46683: 30 27 | 36 3
    955016: 89 63 | 98 24
    515375: 71 54 | 80 15
    314496: 66 30 | 68 4
    65728: 33 31 | 40 12
    443889: 73 38 | 76 17
    704977: 86 41 | 89 2
    327763: 58 51 | 67 30
    4104: 15 9 | 16 2
    171288: 54 24 | 55 17
    195841: 57 22 | 58 9
    262656: 60 36 | 64 8
    165464: 48 38 | 54 20
    134379: 43 38 | 51 12
    994688: 92 60 | 99 29
    684019: 75 64 | 82 51
    593047: 70 63 | 84 7
    920673: 96 33 | 97 20
    885248: 80 72 | 96 8
    373464: 60 54 | 72 6
    439101: 69 48 | 76 5
    402597: 61 56 | 69 42
    110656: 40 36 | 48 4
    110808: 45 27 | 48 6
    513000: 75 45 | 80 10
    842751: 84 63 | 94 23
    984067: 92 59 | 98 35
    886464: 90 54 | 96 12
    1009736: 93 59 | 96 50
    13832: 20 18 | 24 2
    525824: 66 62 | 80 24
    1016496: 90 66 | 97 47
    32832: 30 18 | 32 4
    805688: 92 30 | 93 11
    20683: 24 19 | 27 10
    513856: 72 52 | 78 34
    64232: 36 26 | 39 17
    40033: 33 16 | 34 9
    1729: 10 9 | 12 1
    39312: 33 15 | 34 2
    149389: 50 29 | 53 8
    320264: 66 32 | 68 18
    216027: 59 22 | 60 3

The output is unsorted, as a Tcl array is not sorted, but that is easily fixed, if that is needed. The most important thing is that we can use an array to collect things together without any effort, resulting in a very concise algorithm.

The limit of 100 is not large enough to produce the smallest number that be written as the sum of positive cuvbes in three ways: 87539319 (cf [L2 ])

Functional composition

In real and complex analysis you often use composition of functions, for instance, the function f, defined as:

    f(x) = exp(x^2)

is a composition of the exponential function and a quadratic polynomial.

The [expr] command supports this just like most programming languages:

    proc f {x} {
        return [expr {exp($x*$x)}]
    }

Or, using the new exponentiation operator introduced in Tcl 8.5:

    proc f {x} {
        return [expr {exp($x**2)}]
    }

Sometimes, however, you need a more complicated or more general mechanism. We then need to construct a new function h from two given functions f and g:

    g(x) = f(g(x))

That is very simple to arrange in Tcl:

    proc o {f g} {
        global count

        #
        # Create a new function name
        if { ![info exist count] } {
            set count 0
        } else {
            incr count
        }

        set name "__new__$count"

        proc $name {x} [string map [list F $f G $g] {return [F [G $x]]}]

        return $name
    }

Thus, with:

    proc f {x} {
        return [expr {exp($x)}]
    }
    proc g {x} {
        return [expr {$x*$x}]
    }

we can now create the function "f o g", as:

    % set h [o f g]
    % puts [$h 3]
    ==> 8103.08392758

The function/procedure we just created will work as long as the functions f and g exist. If the new function should be completely independent of these two functions, then we could use the [info body] and other commands to retrieve the definition of these two functions and construct a new function that is independent of that.

Numerical differentiation

The technique of constructing a new procedure can be used in a different context as well. Suppose we need the derivative of a function. While it is perfectly possible to do this manually, it is also error-prone and tedious. Differentiating a function like:

            1 + x^2 + sin(ax)
    f(x) =  -----------------
                 x+ln(x)

is likely to fail the first time.

But we can get a decent numerical approximation without knowing anything about the function:

                    f(x+dx)-f(x)
    f'(x) =   lim   ------------
             dx->0       dx

All we have to do is determine the right value of dx. The literature on numerical analysis suggests to use a step proportional to sqrt(eps), where eps is the smallest number such that 1+eps is distinguishable from 1 (using the limited precision inherent in computer systems).

Of course the size of x plays a role too: x+dx must be different than x in its turn. But for values of x in the "ordinary" range, say 0 to 1000 on a linear scale, where values like 0.0001 are almost zero, sqrt(eps) will do just fine.

eps is roughly 1.0e-12 for most computer systems with double precision floating-point numbers. So:

    proc diff {function} {
        proc $function' {x} [string map [list F $function] {
            set x2 [expr {$x+1.0e-6}]
            return [expr {([F $x2]-[F $x])/1.0e-6}]
        }]
   }

is a procedure that, given a function f, constructs a function f' that returns the (numerical approximation of the) derivative of f. (As there is a natural name for this function, we do not bother to use a scheme like the previous one.)

If you run the following commands:

   proc f {x} {expr {$x*$x}}
   diff f
   puts "1 [f 1] - [f' 1]"
   puts "2 [f 2] - [f' 2]"
   puts "200 [f 200] - [f' 200]"

the result is:

   1 1 - 2.00000099992
   2 4 - 4.00000100065
   200 400 - 399.999997171

so the result is not exact, but will do for many practical purposes.

A little geometry

One fascinating area of mathematics is the construction of all manner of curves: there are many methods to construct one curve from another one and the result is something surprising. Take the family of curves that are formed by circles that are rolled over other circles - the curves you can draw with a "spirograph". The simplest such curve is the cycloid, which is constructed as the trajectory of a fixed point on the perimeter of a circle that is rolled along a straight line.

The script below illustrates this by drawing a seqence of frames from this process:

    proc drawCycloid {cnv radius} {


        set width  [$cnv cget -width]
        set height [$cnv cget -height]

        set steps [expr {($width-1.5*$radius)/10.0}]  ;# Steps of ten pixels

        set xpos0 [expr {1.5*$radius}]
        set xdelt 10.0
        set pi    3.1415926

        #
        # The list of coordinates
        #
        set coords {}

        for { set i 0 } { $i < $steps } { incr i } {

            $cnv delete all

            set xcentr [expr {$xpos0  + $i * $xdelt}]
            set ycentr [expr {$height - $radius}]
            set thetap [expr {1.5*$pi-($i*$xdelt)/$radius}] ;# Angle for the point on the circle
                                                            ;# Clockwise motion

            set xp     [expr {$xcentr         + $radius * cos($thetap)}]
            set yp     [expr {$height-$radius - $radius * sin($thetap)}] ;# Reverse y-coordinate

            lappend coords $xp $yp

            #
            # Draw the circle
            #
            set xcirc1 [expr {$xcentr-$radius}]
            set xcirc2 [expr {$xcentr+$radius}]
            set ycirc1 [expr {$height-2.0*$radius}]
            set ycirc2 $height

            $cnv create oval $xcirc1 $ycirc1 $xcirc2 $ycirc2

            #
            # Draw the spoke
            #

            $cnv create line $xcentr $ycentr $xp $yp

            #
            # Draw the cycloid
            #
            if { $i > 1 } {
                $cnv create line $coords -fill red -width 2
            }

            #
            # Pause ...
            #
            after 100 {set ::next 1}
            vwait ::next
        }
    }

    pack [canvas .c -width 800 -height 300]
    drawCycloid .c 50

Computing the coordinates is straightforward, but you have to remember that:

The y-coordinate on the canvas is upside-down (the origin is in the top-left corner). This is a feature found in many toolkits - for lines of text it does make sense.
The point is rotating in a clockwise fashion, as the circle is rolled to the right.

By pausing 100 ms (the [after] command together with [vwait]) we create the illusion of animation. The [vwait] command makes sure that during the time of waiting, events are processed, even though the execution of the procedure waits until the after command "fires"

Recursive functions

Many interesting series of numbers or families of functions can be defined via recursive relations. For instance the Legendre polynomials can be defined via the recursive relation:

    (n+1) P   (x) - (2n+1) x P (x) + n P   (x) = 0
           n+1                n         n-1

    P (x) = 1, P (x) = x
     0          1

Another famous example, which we will examine in more detail, is the Fibonacci series:

    F(n+2) = F(n+1) + F(n)

    F(0) = 1, F(1) = 1

We can compute these numbers using a recursive procedure:

    proc Fib {n} {

        if { $n < 2 } {
            return 1
        } else {
            set nm2 [expr {$n-2}]
            set nm1 [expr {$n-1}]
            return [expr {[Fib $nm1] + [Fib $nm2]}]
        }
    }

Simple and effective.

Unfortunately, we are doing a lot of unnecessary work, as you can see from the output of a slightly modified version (add: puts "Fib: $n" to see what the argument is):

    % Fib 5
    Fib: 5
    Fib: 4
    Fib: 3
    Fib: 2
    Fib: 1
    Fib: 0
    Fib: 1
    Fib: 2
    Fib: 1
    Fib: 0
    Fib: 3
    Fib: 2
    Fib: 1
    Fib: 0
    Fib: 1
    8

We compute F(0) three times, F(1) five times, F(2) three times, F(3) two times and F(4) one time. This is clearly not very efficient, especially with larger values of n.

We can solve this in two ways:

Use an iterative method instead of a recursive procedure
Use a technique called memoization

The iterative solution is, in this case, fairly simple:

    proc Fib {n} {
        if { $n < 2 } {
            return 1
        } else {
            set Fnm2 1
            set Fnm1 1
            for { set i 2 } { $i <= $n } { incr i } {
                set Fn [expr {$Fnm2+$Fnm1}]
                set Fnm2 $Fnm1
                set Fnm1 $Fn
            }
            return $Fn
        }
    }

But memoization is the more interesting technique, especially as we can do this without any change to the basic algorithm:

    proc memoize {procname} {
        set args \$[join [info args $procname] " $"]

        rename $procname _$procname
        proc $procname [info args _$procname] \
            [string map [list BODY [info body _$procname] \
                              ARGS $args                 \
                              PROC $procname]            \
                {
                   global PROC
                   if { ! [info exists PROC(ARGS)] } {
                       set PROC(ARGS) [_PROC ARGS]
                   }

                   return $PROC(ARGS)
                }]
    }

    memoize Fib

What happens is this:

We rename the original procedure Fib to _Fib
Using the [info] command we can get the definition of the Fib procedure and we use that to build a new procedure.
The [string map] command is very useful for replacing fixed strings
This new procedure first checks if we have already computed the value for the given arguments
If so, it returns that value. If not, it calls _Fib to compute it.

Here is the transformed procedure:

    proc Fib {n} {
        global Fib
        if { ! [info exists Fib($n)] } {
        } else {
            set Fib($n) [_Fib $n]
        }

        return [set Fib($n)]
    }

Running the command "Fib 5" now gives:

    % Fib 5
    Fib: 5
    Fib: 4
    Fib: 3
    Fib: 2
    Fib: 1
    Fib: 0

As you can see, Fib is called only once for any argument.

The [expr] command

The one command for numerical computations in Tcl is [expr]. It supports the usual mathematical operations and functions on both integer and floating-point numbers. It is very much like C or Fortran in that respect, including all the caveats you need to watch out for when dealing with floating-point numbers.

But Tcl has a few extra quirks, you might say, due to the "polymorphic" nature of data in Tcl. Consider a simple procedure to compute the root of a linear equation:

    proc rootLin {a b} {
        #
        # Compute the root x of: a*x + b = 0
        #

        return [expr {-$b/$a}]
    }

This is not a very sophisticated exercise, of course, but you may be surprised to see the result of these commands:

    puts "Root of 2.4x + 3 = 0 is [rootLin 2.4 3]"
    puts "Root of 3x + 2.4 = 0 is [rootLin 3 2.4]"
    puts "Root of 3x + 2   = 0 is [rootLin 3 2]"

The result is:

    Root of 2.4x + 3 = 0 is -1.25
    Root of 3x + 2.4 = 0 is -0.8
    Root of 3x + 2   = 0 is -1

The last answer clearly is incorrect! Why? Because the values 3 and 2 are both integer. That results internally in an integer division which can not give the correct answer -0.66666...

The solution is to redefine the procedure:

    proc rootLin {a b} {
        #
        # Compute the root x of: a*x + b = 0
        # Be careful to convert "$a" to a floating-point value
        #

        return [expr {-$b/double($a)}]
    }

The alternative to use [rootLin 3.0 2] is not acceptable, as you make the user responsable for getting it right, and there might be other computations inside the procedure that he/she has no influence over.

References

The Wiki provides a wealth of mathematically inspired pages:

...

Besides these pages, the Tcllib library contains an extensive math module with packages for statistical computations, linear algebra, integration and interpolation, combinatorics and special functions.

See: [L3 ] for more information.

(Note to self: References are still missing)

AM (1 september 2008) Here is a little application to create those nice mathematical animations you can find on mathworld.wolfram.com - Creating Mathematical Animations

Category Tutorial

Category Mathematics