Open ivan-pi opened 2 years ago
Maybe i can just add this:
pure logical function opposite_sign(x1,x2)
real(wp),intent(in) :: x1,x2
opposite_sign = sign(1.0_wp,x1) /= sign(1.0_wp,x2)
end function opposite_sign
The isign
function was from toms748
. I need to check if there is some subtle difference because they are specifically handing the 0 case...
I noticed SciPy actually had the same issue (https://github.com/scipy/scipy/issues/13737) and they introduced a sign-bit comparison to avoid the test failing due to underflow/overflow.
https://github.com/jacobwilliams/roots-fortran/blob/60524bc4f491a50ae1b2e538e1e486551812b9db/src/root_module.F90#L1837
I noticed the sign comparison is performed inconsistently; i.e. in some places as
fa*fc < 0
and in some places you extract the sign explicitly like above.The latter option is in agreement with what Kahaner, Moler, & Nash have to say about this:
The original SLATEC
FZERO
from Shampine & Watts does it like this:Leaving the underflow issues aside, the
sign
version seems to generate a few instructions less than your branchingisign
function. Moreover,gfortran
andifx
produce branchless instructions when using thesign
intrinsic: https://godbolt.org/z/E41P9c4MzAs a matter of clarity, the sign comparison could also be abstracted as a function or operator, e.g.
or
but I guess this is more a matter of taste (I'm not advocating for it, but the compilers seem able to inline such operators just fine).