math

math is a Tcllib module.

Documentation

math::bigfloat
math::bignum
math::calculus
math::calculus
math::combinatorics
math::complexnumbers
math::constants
math::fourier
math::fuzzy
math::geometry
math::interpolate
math::linearalgebra
math::optimize
math::polynomials
math::rationalfunctions
math::roman
math::romberg
math::special
math::statistics

Contents

::math::cov
coefficient of variation of 3 or more arguments
::math::fibonacci
Return the nth Fibonacci number
::math::integrate
calculate the area under a curve defined by an argument of a list of 5 or more x,y data pairs
::math::max
Return the maximum of two or more arguments
::math::mean
Return the arithmetic mean of two or more arguments
::math::min
Return the minimum of two or more arguments
::math::product
Return the product of two or more arguments
::math::random
A simple wrapper over the rand() function of expr, and subject to overflow. Return a random number, value chosen based on an argument selecting one of these ranges:
              (null)          choose a number between 0 and 1
              val             choose a number between 0 and val
              val1 val2       choose a number between val1 and val2
::math::roman
Handling of Roman Numerals.
::math::sigma
population standard deviation value (from 3 or more arguments)
::math::stats
calculate mean, standard deviation, and coefficient of variation from 3 or more arguments.
::math::sum
calculate arithmetic sum of two or more arguments

If you pick tcllib up from its repository, you will be able to get access to even more functionality:

  • ::math::bigfloat::abs, acos, add, addInt2Float, asin, atan, ceil, ...
  • ::math::bignum::cmp, tostr, add, div, fromstr, max, min, mul, abs, iszero, setsign, ...
  • ::math::calculus::eulerStep, heunStep, rungeKuttaStep, ...
  • ::math::combinatorics Beta, factorialList, pascal, InitializeFactorial, InitializePascal, choose, ln_Gamma, factorial
  • ::math::complexnumbers::complex, +, -, *, / mod, pow, real, sin, sqrt, ...
  • ::math::constants::constants::pi e, ln10 onethird eps print-constants, find_huge, find_tiny
  • ::math::fourier::dft, ...
  • ::math::geometry::calculateDistanceToLine, findClosestPointOnLine, angle, areaPolygon, bbox, calculateDistanceToLine, ...
  • ::math::interpolate::interp-linear, neville, ...
  • ::math::linearalgebra::angle, ...
  • ::math::cov, expectDouble, expectInteger, fibonacci, integrate, max, mean, min, product, random, sigma, stats, sum, ...
  • ::math::optimize::minimum, maximum, min_bound_1d, ...
  • ::math::polynomials::polynomial, addPolyn, subPolyn, ...
  • ::math::rationalnumbers::rationalFunction, ratioCmd, evalRatio, ...
  • ::math::special::elliptic_K, ...
  • ::math::special::exponential_Ei, ...
  • ::math::special::Gamma, ...
  • ::math::special::I_n, J-1/2, J0, J1, J1/2, Jn ...
  • ::math::special::legendre, ...
  • ::math::statistics::BasicStats, histogram, corr, ...
  • ::math::statistics::filter, map, samplescount, ...
  • ::math::statistics::Inverse-cdf-exponential, cdf-normal, ...
  • ::math::statistics::pdf-normal, pdf-uniform, ...
  • ::math::statistics::plot-scale, plot-xydata, ...
  • ::math::statistics::random-exponential, random-normal, ...

TR 2010-09-20: I was missing the Kruskal-Wallis test in math::statistics, and since I needed this test, I went to implementing it. I have posted a corresponding feature request on sf: [L1 ]


FF 2007-06-11: I implemented some missing features of math::linearalgebra. (If they are going to be included in tcllib -I hope for that- you can delete this comment.)

LV 2007-06-11 To get code included in tcllib, visit http://tcllib.sf.net/ and submit a feature request from the sf.net project page.


package require math::linearalgebra

namespace eval ::math::linearalgebra {
    namespace export joinMatrix cancelrow cancelcol
    namespace export cofactor cofactorMatrix
    namespace export adjointMatrix
    namespace export determinant invert rank
    namespace export mkRandom
    namespace export round
}


# joinMatrix --
#     Join two matrices along the specified dimension
# Arguments:
#     dim        Dimension to join [1,2]
#     mv1,mv2    Matrices to be treated
#
# Result:
#     Matrix result of the join
#
proc ::math::linearalgebra::joinMatrix {dim mv1 mv2} {
    set result $mv1
    switch -exact -- $dim {
        1 {
            foreach r $mv2 {lappend $result $r}
        }
        2 {
            set l [llength $mv1]
            for {set i 0} {$i < $l} {incr i} {
                lset result $i [concat [lindex $mv1 $i] [lindex $mv2 $i]]
            }
        }
        default {
            return -code error "do we work on ${dim}D matrices?"
        }
    }
    return $result
}

# cancelrow --
#     Return a matrix minus the specified row
# Arguments:
#     matrix     Matrix to work on
#     row        to delete
#
# Result:
#     Matrix result of the operation
#
proc ::math::linearalgebra::cancelrow {matrix row} {
    return [concat [lrange $matrix 0 $row-1] [lrange $matrix $row+1 end]]
}

# cancelrow --
#     Return a matrix minus the specified column
# Arguments:
#     matrix     Matrix to work on
#     col        to delete
#
# Result:
#     Matrix result of the operation
#
proc ::math::linearalgebra::cancelcol {matrix col} {
    set result {}
    foreach r $matrix {
        lappend result [concat [lrange $r 0 $col-1] [lrange $r $col+1 end]]
    }
    return $result
}

# cofactor --
#     Compute the cofactor of the specified element
# Arguments:
#     matrix     Matrix to work on
#     row        
#     col        
#
# Result:
#     The cofactor of element at a{row,col}
#
proc ::math::linearalgebra::cofactor {matrix row col} {
    set submatrix [cancelcol [cancelrow $matrix $row] $col]
    return [determinant $submatrix]
}

# cofactorMatrix --
#     Compute the cofactor of the specified matrix
# Arguments:
#     matrix     Matrix to work on
#
# Result:
#     The cofactor of specified matrix
#
proc ::math::linearalgebra::cofactorMatrix {matrix} {
    set result {}
    lassign [shape $matrix] rows cols
    for {set i 0} {$i < $rows} {incr i} {
        set newrow {}
        for {set j 0} {$j < $cols} {incr j} {
            lappend newrow [expr ((($i+$j)%2)==0?1:-1)*[cofactor $matrix $i $j]]
        }
        lappend result $newrow
    }
    return $result
}

# adjointMatrix --
#     The adjoint matrix is the transpose of the cofactor matrix
# Arguments:
#     matrix
#
# Result:
#     The adjoint matrix
#
proc ::math::linearalgebra::adjointMatrix matrix {
    return [transpose [cofactorMatrix $matrix]]
}

# determinant --
#     Compute the cofactor of the specified element
# Arguments:
#     matrix     Matrix to work on
#     row        
#     col        
#
# Result:
#     The cofactor of element at a{row,col}
#
proc ::math::linearalgebra::determinant matrix {
    set shape [shape $matrix]
    if {[lindex $shape 0] != [lindex $shape 1]} {
        return -code error "determinant only defined for a square matrix"
    }

    switch -exact -- [join $shape x] {
        1x1 {
            return [lindex [lindex $matrix 0] 0]
        }
        2x2 {
            lassign [getrow $matrix 0] a b
            lassign [getrow $matrix 1] c d
            return [expr {($a*$d)-($c*$b)}]
        }
        default {
            set det 0
            set sign 0
            set row_no 0
            foreach row $matrix {
                set sign [expr {($sign==1)?(-1):(1)}]
                set det [expr {
                    $det+$sign*[lindex $row 0]*[cofactor $matrix $row_no 0]
                }]
                incr row_no
            }
            return $det
        }
    }
}

# invert --
#     Perform the matrix inversion using the adjoint method
#     Note: probably the Gauss-Jordan elimination over an
#           augmented matrix is *much* faster.
# Arguments:
#     matrix     Matrix to work on
#
# Result:
#     The inverted matrix
#
proc ::math::linearalgebra::invert matrix {
    set shape [shape $matrix]

    set det [determinant $matrix]
    if {$det == 0} {
        return -code error "cannot invert a singular matrix"
    }

    switch -exact -- [join $shape x] {
        2x2 {
           set temp [mkMatrix 2 2]
           setelem temp 0 0 [getelem $matrix 1 1]
           setelem temp 1 1 [getelem $matrix 0 0]
           setelem temp 0 1 [expr {-1*[getelem $matrix 0 1]}]
           setelem temp 1 0 [expr {-1*[getelem $matrix 1 0]}]
        }
        default {
           set temp [adjointMatrix $matrix]
        }
    }

    return [scale_mat [expr 1./$det] $temp]
}

Useful for testing:

proc ::math::linearalgebra::mkRandom { rows {cols {}} args} {
    if {$cols eq {}} {set cols $rows}
    set min 0
    set max 1
    if {[llength $args] == 1} {
        set max $args
    } elseif {[llength $args] == 2} {
        lassign $args min max
    }
    set result {}
    for {set i 0} {$i < $rows} {incr i} {
       set row {}
       for {set j 0} {$j < $cols} {incr j} {
           lappend row [expr {$min+$max*rand()}]
       }
       lappend result $row
    }
    return $result
}

proc ::math::linearalgebra::round matrix {
    set result {}
    foreach row $matrix {
        set newrow {}
        foreach el $row {
            lappend newrow [expr {round($el)}]
        }
        lappend result $newrow
    }
    return $result
}

Also I believe these should go in math::combinatorics:

proc ::math::combinations {n k} {
    set l {}
    for {set i 0} {$i < $n} {incr i} {lappend l $i}
    return [lcombinations $l $k]
}

proc ::math::lcombinations {list size} {
    if {$size == 0} {return {{}}}
    set retval {}
    for {set i 0} {($i + $size) <= [llength $list]} {incr i} {
        set firstElement [lindex $list $i]
        set remainingElements [lrange $list [expr {$i+1}] end]
        foreach subset [lcombinations $remainingElements [expr {$size-1}]] {
            lappend retval [linsert $subset 0 $firstElement]
        }
    }
    return $retval
}

Here's a possible implementation of rank:

proc ::math::linearalgebra::rank {matrix {checkR {}}} {
    lassign [shape $matrix] rows cols
    if {$checkR eq {}} {set checkR [expr ($cols>$rows)?$rows:$cols]}
    if {$checkR == 1} {return 1}
    set col_idx [combinations $cols $checkR]
    set row_idx [combinations $rows $checkR]
    foreach coll $col_idx {
        foreach rowl $row_idx {
            set m {}
            for {set i 0} {$i < [llength $rowl]} {incr i} {
                set row {}
                for {set j 0} {$j < [llength $coll]} {incr j} {
                    lappend row [getelem $matrix [lindex $rowl $i] [lindex $coll $j]]
                }
                lappend m $row
            }
            set det [determinant $m]
            if {$det != 0} {return $checkR}
        }
    }
    return [rank $matrix [incr checkR -1]]
}

Lars H: That's horribly inefficient (exponential time), unfortunately. The normal way of computing the rank of a matrix is rather to do some kind of factorisation that reveals the rank (e.g. LU factorisation = Gauss elimination, QR factorisation, etc.). With the available facilities, I would suggest doing a SVD decomposition and count nonzero elements of the S vector.

See Also

Bernoulli using math::calculus
Performance and simplification of code - implementing the basic statistics procedure
matrix , by Jonas Beskow
vector , by Jonas Beskow
tcllib/struct/matrix
NAP
VKIT
BLT vector
Additional math functions
the corresponding Wiki sandbox.