# Outer product of tensors



Eval_outer = ->
  p1 = cdr(p1)
  push(car(p1))
  Eval()
  p1 = cdr(p1)
  while (iscons(p1))
    push(car(p1))
    Eval()
    outer()
    p1 = cdr(p1)

outer = ->
  save()
  p2 = pop()
  p1 = pop()
  if (istensor(p1) && istensor(p2))
    yyouter()
  else
    push(p1)
    push(p2)
    if (istensor(p1))
      tensor_times_scalar()
    else if (istensor(p2))
      scalar_times_tensor()
    else
      multiply()
  restore()

yyouter = ->
  i = 0
  j = 0
  k = 0
  ndim = 0
  nelem = 0

  ndim = p1.tensor.ndim + p2.tensor.ndim

  if (ndim > MAXDIM)
    stop("outer: rank of result exceeds maximum")

  nelem = p1.tensor.nelem * p2.tensor.nelem

  p3 = alloc_tensor(nelem)

  p3.tensor.ndim = ndim

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

  j = i

  for i in [0...p2.tensor.ndim]
    p3.tensor.dim[j + i] = p2.tensor.dim[i]

  k = 0

  for i in [0...p1.tensor.nelem]
    for j in [0...p2.tensor.nelem]
      push(p1.tensor.elem[i])
      push(p2.tensor.elem[j])
      multiply()
      p3.tensor.elem[k++] = pop()

  push(p3)