# Evaluate WHT by definition.  Used for the base case of the WHT recursion
# Inputs:
#  n : log of the size of the transform.  N = 2^n.
#  x : input array
#  base : starting index of subvector of x
#  stride : access stride for subvector of x
# Side Effect:
#  The subvector x^N_{base,stride} is modified to be equal to WHT_N x^N_{base,stride}
#  where x^N_{base,stride} = [x(base),x(base+stride),...,x(base+(N-1)*stride)].

EvalSmall := proc(n,x,base,stride)  
  local N,W,xt,i;
    N := 2^n;  W := WHT(n);
    xt := linalg[vector](N);
    for i from 0 to N-1 do
      xt[i+1] := x[base+i*stride];
    od;
    xt := evalm(W&*xt);
    for i from 0 to N-1 do
      x[base+i*stride] := xt[i+1];
    od;
end;


# Evaluate WHT.
# Inputs:
#  W : partition tree with root n.
#  x : array
#  base : base index of subarray of x.
#  stride : stride that subarray of x is accessed.
# Side Effect.
#  The subvector x^N_{base,stride} is modified to be equal to WHT_N x^N_{base,stride}
#  where x^N_{base,stride} = [x(base),x(base+stride),...,x(base+(N-1)*stride
EvalWHT := proc(W,x,base,stride)
  local n, N, children, W1, N1, n1, t, R, S, i, j, k;
  n := W[1];  N := 2^n;  children := W[2];  t := nops(children);
  if t = 0 then
    EvalSmall(n,x,base,stride);  return;
  fi;
  R := N;  S := 1;
  for i from t by -1 to 1 do
    W1 := children[i];  n1 := W1[1];  N1 := 2^n1;  R := R/N1;
    for j from 0 to R-1 do
      for k from 0 to S-1 do
        EvalWHT(W1,x,base+(j*N1*S+k)*stride,S*stride);
      od;
    od;
    S := S*N1;
  od;
end;

# Versions of EvalSmall and EvalWHT that produce address traces rather than
# performing the computation.  The sequence of indices that are accessed in
# the input array (called x in EvalSmall and EvalWHT) are returned.

TraceSmall := proc(n,base,stride)  
  local N,i,addr;
  N := 2^n;  
  addr := seq(base+i*stride,i=0..N-1);
  return [addr];
end;


TraceWHT := proc(W,base,stride)  
  local n, N, children, W1, N1, n1, t, R, S, i, j, k, addr;
  n := W[1];  N := 2^n;  children := W[2];  t := nops(children);
  if t = 0 then
    return TraceSmall(n,base,stride); 
  fi;
  addr := [];  R := N;  S := 1;
  for i from t by -1 to 1 do
    W1 := children[i];  n1 := W1[1];  N1 := 2^n1;  R := R/N1;
    for j from 0 to R-1 do
      for k from 0 to S-1 do
        addr := [op(addr),op(TraceWHT(W1,base+(j*N1*S+k)*stride,S*stride))];
      od;
    od;
    S := S*N1;
  od;
  return addr;
end;