# examples
 package require vexpr 1.1
 namespace import vexpr::*
 set v1 { 0 0 1 } 
 set v2 { 0 1 0 } 
 puts stdout "v1 = $v1"
 puts stdout "v2 = $v2"
 set v3 [ vexpr $v1 + $v2 ]
 puts stdout "v3 = v1 + v2 = $v3 "
 puts stdout "v1 + v2 = [ vexpr $v1 + $v2 ]"
 puts stdout "v1 - v2 = [ vexpr $v1 - $v2 ]"
 puts stdout "v1 . v2 = [ vexpr $v1 . $v2 ]"
 puts stdout "v1 X v2 = [ vexpr $v1 X $v2 ]"
 puts stdout "v3 size  =[ vexpr $v3 size ]"
 puts stdout "v3 normalise = [ vexpr $v3 normalise ]"
 puts stdout "3 * v3 =  [ vexpr 3 * $v3 ]"
 puts stdout "v3 / 3.0 =  [ vexpr $v3 / 3.0 ]"
 puts stdout "v3 shift 3 = [ vexpr $v3 shift 3 ]"
 #
 set v3 [ vexpr ( $v1 + $v2 ) ]
 puts stderr "v3 = ($v1 + $v2 ) = $v3"
 set v3 [ vexpr ( 2 * $v1 + $v2 ) ]
 puts stderr "v3 = (2 * $v1 + $v2 ) = $v3"
 set v3 [ vexpr  2 * ( $v1 + $v2 ) ]
 puts stderr "v3 = 2 * ( $v1 + $v2 ) = $v3"
 set v3 [ vexpr  2 * ( $v1 + $v2 ) + $v2 ]
 puts stderr "v3 = 2 * ( $v1 + $v2 ) + $v2 = $v3"
 set v3 [ vexpr  2 * ( $v1 + $v2 ) + ( 3 * $v2 / 3 ) ]
 puts stderr "v3 = 2 * ( $v1 + $v2 ) + ( 3 * $v2 / 3 ) = $v3"
 set v4 [ vexpr 2 * ( $v1 + $v2 ) - $v3 + ( 2 * ( $v1 X $v2 ) - $v2 ) ]
 puts stdout "v4 = 2 * ( $v1 + $v2 ) - $v3 + ( 2 * ( $v1 X $v2 ) - $v2 )= $v4" 

 # vexpr.tcl
 package provide vexpr 1.1
 namespace eval ::vexpr:: {
     namespace export vexpr
     variable argc 
     variable sum
 }
 proc ::vexpr::vexpr { args } {
 #
 # vexpr
 # an expression evaluator for vectors.
 # brackets are used to enforce the order of operation
 # expression is evaluated from left to right
 # operators allowed are;    ( v1 and v2 are vectors, s is a scalar )
 # v1 +  v2     add
 # v1 -  v2     subtract
 # v1 shift s   shift
 # v1 X  v2     cross product
 # v1 .  v2     dot product
 # s  *  v1     scale of v1   ( note order of scalar and vector is important)
 # v1 /  s      scaling of v1
 # v1 size      norm of v1
 # v1 normalise return a normalised v1
 #
  variable argc
  variable sum
  set depth 0               ;# depth is the depth of brackets encountered
  set argc(0)  0            ;# argc holds the number of arguments found so far
  set sum(0,0) ""           ;# sum holds the 3 components of a sum
  set sum(0,1) ""           ;#             a op b
  set sum(0,2) ""           ;#
  foreach element $args {
     if { $element == "(" } { 
        # go to the next level of brackets and initialise
        incr depth
        set argc($depth)  0
        set sum($depth,0) ""
        set sum($depth,1) ""
        set sum($depth,2) ""
     } elseif { $element == ")" } {
        # finished with the present level, evaluate current expression and 
        # move up in depth.  Store the result.  
        set result [ _getResult $depth ]
        incr depth -1
        set sum($depth,$argc($depth)) $result
        _addToList $depth
     } else {
        # store the next argument in the list
        set sum($depth,$argc($depth)) $element
        _addToList $depth
     }
        
  }
  return [ _getResult $depth ]
 }

 proc ::vexpr::_addToList { depth } {
  # 
  # Increment the number of arguments at this depth
  # if we have 3 then evaluate and store
  #
  variable argc
  variable sum
  incr argc($depth)
  if { $argc($depth) > 2 } {
    set sum($depth,0) [ _voperate $sum($depth,0) $sum($depth,1) $sum($depth,2) ]
    set argc($depth) 1
    set sum($depth,1) ""
    set sum($depth,2) ""
  }
 }

 proc ::vexpr::_getResult { depth } {
  # 
  # Return the result.  If we have had 1 binary op
  # it will have been completed, but if unary it will not have been
  #
  variable argc
  variable sum
  if { $sum($depth,1) == "" } {
      set result $sum($depth,0)
  } else {
      set result [ _voperate $sum($depth,0) $sum($depth,1) $sum($depth,2) ]
  }
 }

 proc ::vexpr::_voperate { a op { b {} } } {
  switch -- $op {
  +         { set result [ _vadd    $a $b ] }
  -         { set result [ _vsub    $a $b ] }
  shift     { set result [ _vshift  $a $b ] }
  X         { set result [ _vcross  $a $b ] }
  .         { set result [ _vdot    $a $b ] }
  *         { set result [ _vscale  $a $b ] }
  /         { set result [ _vdivide $a $b ] }
  normalise { set result [ _vnormalise $a ] }
  size      { set result [ _vsize $a ] }
  default   { puts stderr "Unkown operator $op" }
  }
 }

 proc ::vexpr::_vshift { v1 s } {
 #
 # Shift vector by s
 # 
    set r ""
    foreach a $v1  {
     lappend r [ expr { $a + $s } ]
    }
    return $r
 }

 proc ::vexpr::_vadd { v1 v2 } {
 #
 # Add two Vectors
 #
    set r ""
    foreach a $v1  b $v2 {
     lappend r [ expr { $a + $b } ]
    }
    return $r
 }

 proc ::vexpr::_vsub { v1 v2 } {
 #
 # Subtract two Vectors
 #
    set r ""
    foreach a $v1 b $v2 {
     lappend r [ expr { $a - $b } ] 
    }
    return $r
 }

 proc ::vexpr::_vdot { v1 v2 } {
 # 
 # dot product of two Vectors
 #
    set sum 0.0 
    foreach a $v1 b $v2 {
     set sum [ expr { $sum + ( $a * $b ) } ]
    }
    return $sum
 }

 proc ::vexpr::_vnormalise { v1 } {
 #
 # Normalise a vector
 #
    set sum [ _vsize $v1 ]
    return  [ _vdivide $v1 $sum ]
 }

 proc ::vexpr::_vsize { v1 } {
 #
 # size of a vector
 #
    set sum  [ _vdot $v1 $v1 ]
    set sum  [ expr { sqrt ( $sum ) } ]
    return $sum
 }

 proc ::vexpr::_vcross { v1 v2 } {
 #
 # cross product of two Vectors
 #
    set r ""
    set ax [ lindex $v1 0 ]
    set ay [ lindex $v1 1 ]
    set az [ lindex $v1 2 ]
    set bx [ lindex $v2 0 ]
    set by [ lindex $v2 1 ]
    set bz [ lindex $v2 2 ]
    set cx [ expr { ($ay*$bz) - ($az*$by) } ]
    set cy [ expr { ($az*$bx) - ($ax*$bz) } ]
    set cz [ expr { ($ax*$by) - ($ay*$bx) } ]
    set r [ list $cx $cy $cz ]
    return $r
 }

 proc ::vexpr::_vscale { s v1 } {
 #
 # scale a vector 
 #
    set r ""
    foreach a $v1 {
     lappend r [ expr { $a * $s } ]
    }
    return $r
 }

 proc ::vexpr::_vdivide { v1 s } {
 #
 # scale a vector 
 #
    set r ""
    foreach a $v1 {
     lappend r [ expr { $a / $s } ]
    }
    return $r
 }

 foreach {ax ay az} v1 {bx by bz} v2 {break;};

 set ax [lindex $v1 0];
 set ay [lindex $v1 1];
 set az [lindex $v1 2];
 ...

 return [list \
     [expr {($ay*$bz) - ($az*$by)}] \
     [expr {($az*$bx) - ($ax*$bz)}] \
     [expr {($ax*$by) - ($ay*$bx)}] \
 ];

 set cx [ expr { ($ay*$bz) - ($az*$by) } ]
 set cy [ expr { ($az*$bx) - ($ax*$bz) } ]
 set cz [ expr { ($ax*$by) - ($ay*$bx) } ]
 set r [ list $cx $cy $cz ]
 return $r