for i=1:500
  a = rand (2,1);
  b = rand (2,1);
  a1 = norm (a,2);
  b1 = norm (b,2);
%  arg = dot (a,b)/(a1*b1);
  theta = acos (dot (a,b)/(a1*b1));
%  theta = asin (sqrt(1 - (dot (a,b)/(a1*b1))^2));
%  theta = 0.5 * acos (2*(dot (a,b)/(a1*b1))^2 - 1);
  if (abs (subspace (a, b) - theta) > 100*eps)
%  if (arg < 0.999 && abs (subspace (a, b) - theta) > 100*eps)
    abs (subspace (a, b) - theta)
    break;
  end
end

