Fortran - Almost equal

Author

Almost equal

rih5342

2005-09-16, 6:59 pm

What's the best way to do almost equal?

The following code works but I'm wondering if there is a better way. I'm
working with a geometry
database and this is done ~30,000 times, over and over.

EQUAL = .FALSE
IF( X .EQ. Y ) THEN
EQUAL = .TRUE.
ELSEIF( X .EQ. 0. ) THEN
IF(ABS(Y) .LE. EPS) EQUAL = .TRUE.
ELSEIF( Y .EQ. 0. ) THEN
IF(ABS(X) .LE. EPS) EQUAL = .TRUE.
ELSEIF( (X .NE. 0.) .AND. (Y .NE. 0.) ) THEN
DDD = X/Y - 1.0
IF( ABS(DDD) .LE. EPS) EQUAL = .TRUE.
ENDIF

Thanks in advance.
Rob

glen herrmannsfeldt

2005-09-16, 6:59 pm

rih5342 wrote:

> What's the best way to do almost equal?

> The following code works but I'm wondering if there is a better way. I'm
> working with a geometry
> database and this is done ~30,000 times, over and over.

> EQUAL = .FALSE
> IF( X .EQ. Y ) THEN
> EQUAL = .TRUE.
> ELSEIF( X .EQ. 0. ) THEN
> IF(ABS(Y) .LE. EPS) EQUAL = .TRUE.

These are doing absolute error, which is the only thing
you can do when one is zero...

> ELSEIF( Y .EQ. 0. ) THEN
> IF(ABS(X) .LE. EPS) EQUAL = .TRUE.
> ELSEIF( (X .NE. 0.) .AND. (Y .NE. 0.) ) THEN
> DDD = X/Y - 1.0

then this one does relative error.
(Though I probably would have done (X-Y)/Y.)

> IF( ABS(DDD) .LE. EPS) EQUAL = .TRUE.
> ENDIF

does your problem expect relative error or absolute error?
I think the absolute error,

IF(ABS(X-Y)<EPS) is fairly commonly used, much simpler,
but it won't work if you need relative error.

You could avoid the divide problem with

IF(ABS(X-Y)<EPS*ABS(Y)) or maybe something like

IF(ABS(X-Y)<EPS*MAX(1.,ABS(Y))

which I think is more fair than your way of using the
same EPS for both relative and absolute error (*). OR

IF((X-Y)**2<EPS**2*(1+Y**2))

which makes a smooth transition between absolute are
relative errors, and avoids the many conditional tests
that yours has.

Note that on most modern processors conditional jumps
are much more expensive than most arithmetic operations.

(*) Yours has a strange discontinuity. Consider EPS=.001,
X=1e-20, Y=1e-40. It will fail your relative error test,
but when Y decreases to 0.0 it will pass the absolute
error test.

-- glen

Paul Van Delst

2005-09-16, 6:59 pm

rih5342 wrote:
> What's the best way to do almost equal?
>
> The following code works but I'm wondering if there is a better way. I'm
> working with a geometry
> database and this is done ~30,000 times, over and over.
>
> EQUAL = .FALSE
> IF( X .EQ. Y ) THEN
> EQUAL = .TRUE.
> ELSEIF( X .EQ. 0. ) THEN
> IF(ABS(Y) .LE. EPS) EQUAL = .TRUE.
> ELSEIF( Y .EQ. 0. ) THEN
> IF(ABS(X) .LE. EPS) EQUAL = .TRUE.
> ELSEIF( (X .NE. 0.) .AND. (Y .NE. 0.) ) THEN
> DDD = X/Y - 1.0
> IF( ABS(DDD) .LE. EPS) EQUAL = .TRUE.
> ENDIF

I use

ABS( x - y ) < ( ULP * SPACING( MAX(ABS(x),ABS(y)) ) )

The explanation from my source code is:

! PROCEDURE:
! The test performed is
!
! ABS( x - y ) < ( ULP * SPACING( MAX(ABS(x),ABS(y)) ) )
!
! If the result is .TRUE., the numbers are considered equal.
!
! The intrinsic function SPACING(x) returns the absolute spacing of numbers
! near the value of x,
!
! { EXPONENT(x)-DIGITS(x)
! { 2.0 for x /= 0
! SPACING(x) = {
! {
! { TINY(x) for x == 0
!
! The ULP optional argument scales the comparison.
!
! James Van Buskirk and James Giles suggested this method for floating
! point comparisons in the comp.lang.fortran newsgroup.

where "ULP" is a user supplied scale factor (for those case where your numbers are close, but not
that close). Definition of ULP:

! OPTIONAL INPUT ARGUMENTS:
! ULP: Unit of data precision. The acronym stands for "unit in
! the last place," the smallest possible increment or decrement
! that can be made using a machine's floating point arithmetic.
! A 0.5 ulp maximum error is the best you could hope for, since
! this corresponds to always rounding to the nearest representable
! floating-point number. Value must be positive - if a negative
! value is supplied, the absolute value is used.
! If not specified, the default value is 1.
! UNITS: N/A
! TYPE: INTEGER
! DIMENSION: Scalar
! ATTRIBUTES: OPTIONAL, INTENT( IN )
!

cheers,

paulv

--
Paul van Delst
CIMSS @ NOAA/NCEP/EMC

Herman D. Knoble

2005-09-19, 7:58 am

Rob: What you need is a way to scale the comparison to the unit
interval and then compare using the machine comparison tolerance,
Epsilon(x). That is:

abs(x-y) <= max(abs(x),abs(y))*fuzz
where fuzz is, for example, 2*Epsilon(x)

This is illustrated with some other "tolerant" or "fuzzy"
functions in the code:
http://ftp.cac.psu.edu/pub/ger/fortran/hdk/eps.f90

Skip Knoble

On Fri, 16 Sep 2005 17:04:46 -0400, "rih5342" <rob.hickey@nospam.alum.mit.edu> wrote:

-|What's the best way to do almost equal?
-|
-|The following code works but I'm wondering if there is a better way. I'm
-|working with a geometry
-|database and this is done ~30,000 times, over and over.
-|
-|EQUAL = .FALSE
-|IF( X .EQ. Y ) THEN
-| EQUAL = .TRUE.
-| ELSEIF( X .EQ. 0. ) THEN
-| IF(ABS(Y) .LE. EPS) EQUAL = .TRUE.
-| ELSEIF( Y .EQ. 0. ) THEN
-| IF(ABS(X) .LE. EPS) EQUAL = .TRUE.
-| ELSEIF( (X .NE. 0.) .AND. (Y .NE. 0.) ) THEN
-| DDD = X/Y - 1.0
-| IF( ABS(DDD) .LE. EPS) EQUAL = .TRUE.
-| ENDIF
-|
-|Thanks in advance.
-|Rob
-|
-|

John H. Lindsay

2005-09-19, 6:58 pm

Hi Skip & other folks,

One function that I have found handy here is my
Relative Difference function RELDIF:

Function RelDif ( x, y )

If ( x .eq. y ) Then
RelDif = 0.0
Else
RelDif = 2.0 * Abs ( x - y ) /
* Abs ( x ) + Abs ( y )
End if

End

Please excuse my FORTRAN IV (or Fortran 77)!

John.

Herman D. Knoble wrote:
> Rob: What you need is a way to scale the comparison to the unit
> interval and then compare using the machine comparison tolerance,
> Epsilon(x). That is:
>
> abs(x-y) <= max(abs(x),abs(y))*fuzz
> where fuzz is, for example, 2*Epsilon(x)
.......

rih5342

2005-09-19, 6:58 pm

Hi Glen, Paul, and Herman.

I guess I should have searched on CLF first.

In any case I now have your useful posts and source suggestions to nail
this down.

Thank you.

Carlie J. Coats

2005-09-20, 7:02 pm

John H. Lindsay wrote:
> Hi Skip & other folks,
>
> One function that I have found handy here is my
> Relative Difference function RELDIF:
>
> Function RelDif ( x, y )
>
> If ( x .eq. y ) Then
> RelDif = 0.0
> Else
> RelDif = 2.0 * Abs ( x - y ) /
> * Abs ( x ) + Abs ( y )
> End if
>
> End
>
> Please excuse my FORTRAN IV (or Fortran 77)!
>
> John.

If you're really interested in performance here, there are a couple
of problems with this formulation on most modern hardware:

absolute value may be slow on hardware that doesn't have
predication (the comparisons and branches don't play friendly
with deeply pipelined CPUs) [If the compilers are good enough,
this objection does not hold for POWER, Itanium, and Alpha.]

Divisions are slow on almost everyone's hardware.

For the environmental modeling work I do (and, pragmatically, for the
tolerances implied by its data), and using the fact that squaring has
the same order properties as absolute value, I tend to use the following
formulation for "certainly not equal":

LOGICAL FLTERR
REAL PP, QQ

FLTERR( PP, QQ ) =
& ( (PP - QQ)**2 .GT. 1.0E-10*( PP*PP + QQ*QQ + 1.0E-20 ) )

The "+1.0E-20" term avoids the zeros problem; the available fine
grained parallelism is maximized; there is no (slow) division, and
only one (pipeline-breaking) comparison...

--

Carlie J. Coats, Jr. carlie.coats@baronams.com
Environmental Modeling Center carlie_coats@ncsu.edu
Baron Advanced Meteorological Systems

Sponsored Links