3

The Netlib library module for the Brent root-finding function checks that the signs of two variables are different like this:

if (fa * (fb/dabs(fb)) .le. 0.0d0) go to 20
...

Why would this check include /dabs(fb) instead of being simply (fa*fb) .le. 0.0d0? I did a quick check with Python and it seems for very large values (+/-1e200) for x and y, where x*y=+/- inf, the comparison x*y<= 0 still works correctly.

iboisver
  • 431
  • 4
  • 13
  • Q: If all you're interested in is the sign ... then there's nothing "wrong" with `(fb/dabs(fb))` ("+1/-1"). I'm guessing it just mitigates the risk of overflow? For historical reasons? – paulsm4 Nov 05 '20 at 22:15
  • @paulsm4 It looks like you might be right about mitigating risk of overflow. [This (old) page](https://docs.oracle.com/cd/E19957-01/805-4940/6j4m1u7pg/index.html) describes default Fortran behaviour of a floating-point overflow (different for f77 and f90) and trapping exceptions. I guess gcc would have similar options. – iboisver Nov 05 '20 at 22:34

1 Answers1

3

Fortran has never specified a function like signs_differ(x,y), so one generally implements such a thing personally.

x*y<0 (and x*y.lt.0) is not asking the same thing as "are x and y of different sign?". While the product of x and y being positive means x and y are the same sign in the (mathematical) real numbers, this is not true for (computational) floating point numbers.

Floating point multiplication x*y may overflow, result in a signed infinite value (raising a IEEE flag) with the comparison returning the expected logical value, but that isn't always true. There were many non-IEEE systems and IEEE systems may see that flag being raised and abort (or have some expensive handling diversion). That's totally not the same thing as "do x and y have the same sign?".

x*(y/dabs(y)) doesn't overflow, is "portable" and is potentially cheaper than (x/dabs(x))*(y/dabs(y)) - ignoring the issues surrounding dabs() and signed zeros.

Modern Fortran has functions such as sign, ieee_copy_sign and ieee_signbit which didn't exist 40 years ago.

francescalus
  • 30,576
  • 16
  • 61
  • 96