In sqrt_Fp2, I think there is a subtle error in the implementation of the square root that leads to wrong answers in the specific case where u = a + 0i, and a is not a quadratic residue modulo p.
The issue arises when calculating (a^2 + b^2)^((p+1)/4). When b = 0, this is a^((p+1)/2) = aa^((p-1)/2)) = aL(a), where L(a) is the Legendre symbol modulo p. This is then added to a on the next line. But when a is a quadratic non-residue, this then sets t_0 to zero, and subsequently therefore the whole square root is set to zero.
It might be that this doesn't occur in practice for some reason but if so I can't see why.
I've written up a more full exposition which contains a (very modest) fix:
squareroots1.pdf
I apologise if this actually isn't an issue and I've missed a mitigation for it.