Floating point equality comparison if(a==0)

From: Dmitri Kavetski (no_at_one.knows)
Date: 03/31/04 

Date: Wed, 31 Mar 2004 21:54:25 +1000

Having looked at a few numerical codes (including netlib, etc.), I noticed
that in many cases they rely on equality-comparing real numbers, eg, in
if(a==const) statements. This is potentially dangerous if the numbers are
put in registers before the comparison (Steve's "perils of real numbers"
springs to mind), but in the single but common case (const=0.0) it is not
clear to me how can this be avoided without considerable extra complexity.

So I would be interested in hearing the opinion of compiler-writers
regarding such comparisons, as well as hearing from numerical programmers
that have written or deal with such codes. Perhaps there are some simple
ways I am unaware of...

If const<>0 then it is easy to replace the equality test with an interval
test with tolerance ~ epsilon(const)*const. However in the special case
const=0.0 the logic fails since any roundoff (even tiny(a)) will negate the
test. The problem is that 0.0 does not have an obvious scale.

So, my question is.... is there any way (or, for that matter, real need) to
fix
!*********
! code 1
a=0.0
!...lots of code Z that may or may not use "a"
if(a==0.0)then
! reliably(?) assume no action on "a" took place
endif
!*********
without introducing additional logical variables such as
!*******
! code 2
a=0.0; aIsZero=.true.
!...lots of code Z, WHICH sets aIsZero=.false. if "a" is used.
if(aIsZero)then
! reliably assume no action on "a" took place
endif
!*******
Note that simply comparing if(abs(a)<=epsilon(a)) is not quite the same,
since it may ignore small but significant changes to a. The problem seems
that the very act of placing the number "a=0" (unlike any other number) in a
higher precision register may perturb it in a way that makes it impossible
(as far as I can tell) to identify whether it was "register junk" or
"genuine computation". Or is there some _guaranteed_ floating point
behaviour that can be relied on (now and in the future) ?

In the above example the penalty for using a logical flag is trivial, but
what if "a" is a large array with many dimensions (indeed tests such as
these can occur frequently in matrix computations). Moreover, what if "a" is
passed through many subroutines - in that case aIsZero has to be dragged
along as well...

** IS THERE a widely-accepted robust and _consistent_ way of implementing
such tests in Fortran? I am also curious to what extent are legacy codes
vulnerable to such perils of more sophisticated floating-point
implementation (with higher-precision registers, etc.) in modern compilers
(eg, Compaq/Intel Fortran). Also, I presume all languages (eg, C++) are
similarly affected.

Dmitri

Relevant Pages

  • Re: Floating point equality comparison if(a==0)
    ... Or is there some _guaranteed_ floating point ... -|implementation (with higher-precision registers, etc.) in modern compilers ... -|(eg, Compaq/Intel Fortran). ... Penn State Information Technology Services ...
    (comp.lang.fortran)
  • Re: builtin lists and Intel SSE support?
    ... However the floating point speedups could also be related to the ... > utilising extra integer registers). ... > ability to intermix integer and floating-point instructions in the same ... > This achieves potentially four times the computational work of x87 ...
    (comp.lang.lisp)
  • Re: performance gain with 32 bit Linux on 64 bit system
    ... I will have a 32 bit Linux running on ... Floating Point has been processed 64bits at a time for a very very long ... added a bunch of additional registers which are only available to a 64 ... Generally if the data set is mostly 32 integers and it's large but not so ...
    (comp.os.linux.misc)
  • Re: Combining two MMX registers into one SSE register?
    ... I already had resigned to the fact not to use MMX under ... Can x87 instructions still be used by 64-bit applications? ... About x87 and mmx registers ... Can floating point registers be used in 64-bit Windows? ...
    (comp.lang.asm.x86)
  • Re: performance gain with 32 bit Linux on 64 bit system
    ... I will have a 32 bit Linux running on ... Floating Point has been processed 64bits at a time for a very very long ... added a bunch of additional registers which are only available to a 64 ... Generally if the data set is mostly 32 integers and it's large but not so ...
    (comp.os.linux.misc)