NANs in precessing HBR final spin fitting formulas

I have seen NANs when using the precessing HBR final spin fit (as implemented in nrutilus.py, using the 'M3J4' variant). This actually comes from trying to debug an alternative implementation of this fit which is used internally in SEOBNRv4P. In case this is user error, I apologise in advance.

In particular, when looking at very large number of cases (10^6), I found O(~few) case where the final spin is NAN. This is coming from the fact that the fit features a square root and for some very particular cases things inside this square root become negative (see Eq(16) in https://arxiv.org/pdf/1605.01938.pdf). Here are 2 concrete examples for which the code produces NANs(all angles are in radians):

m_{1}	m_{2}	\chi_{1}	\chi_{2}	t_{1}	t_{2}	\phi_{12}
63.0361223508	8.42503598697	0.454719654835	0.699876847904	3.08181486287	1.84143980375	3.03597524558
89.9182019495	17.1217349168	0.643333566359	0.47846069899	3.1309081245	1.60082705402	4.39151182036

In both cases, if one takes the absolute value of the argument of the square root one obtains values very close to 0.

Attached is the script that I used to find these cases. minimal_example.py