### contract =====================================================================

Tags
----
scripting, JS, internal, treenode, general concept

Parameters
----------
a,i,j

General description
-------------------
Contract across tensor indices i.e. returns "a" summed over indices i and j.
If i and j are omitted then 1 and 2 are used.
contract(m) is equivalent to the trace of matrix m.

###



Eval_contract = ->
  push(cadr(p1))
  Eval()
  if (cddr(p1) == symbol(NIL))
    push_integer(1)
    push_integer(2)
  else
    push(caddr(p1))
    Eval()
    push(cadddr(p1))
    Eval()
  contract()

contract = ->
  save()
  yycontract()
  restore()

yycontract = ->
  h = 0
  i = 0
  j = 0
  k = 0
  l = 0
  m = 0
  n = 0
  ndim = 0
  nelem = 0
  ai = []
  an = []

  p3 = pop()
  p2 = pop()
  p1 = pop()

  if (!istensor(p1))
    if (!isZeroAtomOrTensor(p1))
      stop("contract: tensor expected, 1st arg is not a tensor")
    push(zero)
    return

  push(p2)
  l = pop_integer()

  push(p3)
  m = pop_integer()

  ndim = p1.tensor.ndim

  if (l < 1 || l > ndim || m < 1 || m > ndim || l == m \
  || p1.tensor.dim[l - 1] != p1.tensor.dim[m - 1])
    stop("contract: index out of range")

  l--
  m--

  n = p1.tensor.dim[l]

  # nelem is the number of elements in "b"

  nelem = 1
  for i in [0...ndim]
    if (i != l && i != m)
      nelem *= p1.tensor.dim[i]

  #console.log "nelem:" + nelem
  p2 = alloc_tensor(nelem)
  #console.log "p2:" + p2

  p2.tensor.ndim = ndim - 2

  j = 0
  for i in [0...ndim]
    if (i != l && i != m)
      p2.tensor.dim[j++] = p1.tensor.dim[i]


  a = p1.tensor.elem
  b = p2.tensor.elem

  #console.log "a: " + a
  #console.log "b: " + b

  for i in [0...ndim]
    ai[i] = 0
    an[i] = p1.tensor.dim[i]

  for i in [0...nelem]
    push(zero)
    for j in [0...n]
      ai[l] = j
      ai[m] = j
      h = 0
      for k in [0...ndim]
        h = (h * an[k]) + ai[k]
      push(a[h])
      #console.log "a[h]: " + a[h]
      add()
      #console.log "tos: " + stack[tos-1]
    b[i] = pop()
    #console.log "b[i]: " + b[i]
    for j in [(ndim - 1)..0]
      if (j == l || j == m)
        continue
      if (++ai[j] < an[j])
        break
      ai[j] = 0

  if (nelem == 1)
    push(b[0])
  else
    push(p2)

  #console.log "returning: " + stack[tos-1]