Robust sign comparison

[ivan-pi]

roots-fortran/src/root_module.F90

Line 1837 in 60524bc

if ((isign(fa)*isign(fc)) < 0) then

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:

IF (SIGN(1.0E0,FB) .NE. SIGN(1.0E0,FC)) ...

Leaving the underflow issues aside, the sign version seems to generate a few instructions less than your branching isign function. Moreover, gfortran and ifx produce branchless instructions when using the sign intrinsic: https://godbolt.org/z/E41P9c4Mz

As a matter of clarity, the sign comparison could also be abstracted as a function or operator, e.g.

if (sne(fb,fc)) ...  ! sign not equal

or

if( fb .oppsign. fc) ...   ! opposite sign

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).

[jacobwilliams]

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...

[ivan-pi]

I noticed SciPy actually had the same issue (scipy/scipy#13737) and they introduced a sign-bit comparison to avoid the test failing due to underflow/overflow.