`calculate_ellipse_parameters` fails to calculate the correct values for `a` and `b` when the ellipse rotation approaches 0.25 or 1.25 about ellipsesampling.jl HOT 1 CLOSED

The culprit is likely to be line 29 in src/matrix_to_parameters.jl which subtracts cos(α)^4 - sin(α)^4. When alpha is 0.25pi and 1.25pi cos(α) is equal to sin(α). However due to numerical issues, both when calculating alpha and when calculating the subtraction of the 4th power of these functions, the actual value is not quite a zero:

println(cos(0.25pi)^4 - sin(0.25pi)^4); println(cos(α)^4 - sin(α)^4)

1.3877787807814457e-16
1.3877787807814457e-16

The behaviour has also been noted for negative values of alpha close to -0.25pi and -1.25pi.

We can fix this issue by enhancing the precision of our calculations using BigFloat on the normalised matrix Hw on line 26. A precision of 64 gives a satisfactory result, however, it does impact the performance of the function by about 40 times. Lower precision values also gave good results for this test case, however, the additional precision was valuable when the magnitude of the two variables is significantly different.

A satisfactory compromise is to selectively enhance the precision if alpha is detected to be close to these problematic values using: isapprox(abs(rem(α/pi, 1)), 0.25, atol=1e-2)

from ellipsesampling.jl.