fortran - 为什么计算两个变量乘积的符号要用绝对值？

Question

Brent root-finding function 的 Netlib 库模块检查两个变量的符号是否不同，如下所示:

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

为什么这个检查会包括 /dabs(fb) 而不是简单的 (fa*fb) .le。 0.0d0？我用 Python 进行了快速检查，似乎 x 和 y 的值非常大 (+/-1e200)，其中 x*y=+/- inf，比较 x*y<= 0 仍然可以正常工作。

Answer 1

Fortran 从未指定过类似signs_differ(x,y) 的函数, 所以一般都是亲自去实现这样的事情。

x*y<0 (和 x*y.lt.0 )与“x 和 y 的符号不同吗？”问的不是同一件事。虽然 x 和 y 的乘积为正意味着 x 和 y 在(数学)实数中符号相同，但对于(计算的) float 而言并非如此。

浮点乘法 x*y可能溢出，导致带符号的无限值(引发 IEEE 标志)，比较返回预期的逻辑值，但这并不总是正确的。有许多非 IEEE 系统和 IEEE 系统可能会看到该标志被升起并中止(或者有一些昂贵的处理转移)。这与“x 和 y 的符号相同吗？”完全不同。

x*(y/dabs(y))不会溢出，“便携”并且可能比 (x/dabs(x))*(y/dabs(y)) 便宜- 忽略围绕 dabs() 的问题和有符号的零。

现代 Fortran 语言有 sign 等函数, ieee_copy_sign和 ieee_signbit这在 40 年前是不存在的。