function [C, pathA, pathB] = lca(A, B, h, G)
%
% [C, pathA, pathB] = lca(A, B, heights, G)
% C is the least common ancestor of A and B in a rooted tree G if
% C is the "lowest possible" ancestor of A and B.
% pathA is the sequence of nodes going from A to C, and similarly for pathB (including endpoints).
% This algorithm takes time linear in |pathA| + |pathB|.
% The height of each node below the root can be computed using dfs.

directed = 0;
assert(acyclic(G,directed)); % ensure G is a tree

pathA = [A];
pathB = [B];

if h(A) < h(B) % B is deeper, so raise it to the level of A
  for i=1:h(B)-h(A)
    ps = parents(G, B);
    assert(length(ps)==1);
    p = ps(1);
    pathB = [pathB p];
    B = p;
  end
elseif h(B) < h(A)
  for i=1:h(A)-h(B)
    ps = parents(G, A);
    assert(length(ps)==1);
    p = ps(1);
    pathA = [pathA p];
    A = p;
  end
end
assert(h(A) == h(B));
if h(A)==0 | A==B
  done = 1;
else
  done = 0;
end
while ~done
  ps = parents(G,A);
  assert(length(ps)==1);
  A = ps(1);
  ps = parents(G,B);
  assert(length(ps)==1);
  B = ps(1);
  if A==B
    done = 1;
  end
  pathA = [pathA A];
  pathB = [pathB B];
end
assert(A==B);
C = A;


