This commit fixes the bug in phs_symmetrization when it fails to produce ph eigenstates
with PH^2 = +1 when one of the wave functions is mapped exactly onto another one by PH,
as in this case the PH symmetrized wf's from the two are identical. This is the case
with wfs = np.eye(2) and  particle_hole = sigma[1], the second vector returned is an
eigenstate with eigenvalue -1 not +1. In general this can be even worse, if there are
multiple such vector pairs some vectors may not be ph eigenstates at all.

The new code uses an improved version of the matrix square root method from
Applied Mathematics and Computation 234 (2014) 380-384.

Tests are modified to clear up the random unitary generation process (using kwant.rmt).
Tests are now guaranteed to feed orthonormal set of vectors to phs_symmetrization().
Add extra tests for the fully off-diagonal case where the old method fails.