def norm(x):
  return x / np.linalg.norm(x)

def procrustes(p1, p2):
  u, _, vh = np.linalg.svd(p2.T@p1)
  if np.linalg.det(u@vh) < 0:
    u = np.stack([u[:, 0], u[:, 1], -u[:, 2]])

  return u@vh

def get_xyz():
  x, y = -1 + 2*np.random.rand(2,3)
  # x, y = norm(x), norm(y)
  z = np.cross(x,y)
  return np.stack([x,y,z], axis=0)

x, y, z = [get_xyz() for _ in range(3)]

# run triplet procrustes
R12 = procrustes(x, y)
R23 = procrustes(y, z)
R31 = procrustes(z, x)

print(R23 @ R12 @ R31)