A = p1'*p2; product = crossdiags(A);