Real number tolerance and if()

[gelarimer]

After reading Doctor Fortran's article on 'The Perils of Real Numbers, Part 1', Iexperimented with some old code and found the following:

REAL Sigma
CHARACTER(50) sStr

if(Sigma .GT. 6.3) RETURN

write(sStr, '(F35.30)') Sigma
..

Program did not RETURNand sStr valuewas
' 6.3000001907348632815', so Sigma was .GT. 6.3.How is it possible for Sigma.GT. 6.3 to be evaluated as .FALSE.?

Thanks for any replies.

[Steven_L_Intel1]

That was Dave Eklund's article, not mine, but..

Perhaps you need to read it again. Decimal fractions are often not exactly representable in binary floating point. 6.3 is one of them. Compile and run this program and then think about your question.

write (*,'(F35.30)') 6.3
end

[greldak]

Also you don't convert until after you have tested sigma. The contents of the character variable are irrelevant - you're not testing that and it has not yet affected the value of the variable you are testing.

What is the value of Sigma when you carry out the test, if its undefined you can expect erratic behaviour, somtimes it will return other times not.

[gelarimer]

Sigma was read fromafilecontaing the text6.3. What I failed to understand was that both the variable Sigma = 6.3and the constant 6.3 are representedthe same in memory and both representationsare not exact, so if(6.3000001907.. .GT. 6.3000001907..) evaluates to .FALSE.

[jimdempseyatthecove]

The number 6.3 cannot be represented using the internal storage format for floating point numbers. This is true regardless of if the number stored is REAL(4), REAL(8), REAL(16). This is due to 0.1 being an infinately long repeting binary fraction. (0.3 will also be infinately long). The internal representation can store a number that is approximately 6.3 to within one bit value of the precision of the storage format (4, 8, 16 byte) for the particular number being represented. This characteristic is similar to rounding but it would more properly be described as uncertanty.

The IF test compares the internal formats of the two numbers (or expressions) under the assumption that there is no uncertanty in the precision. Therefor, depending on how 6.3 is derived in the program you will either end up with 6.3+uncertanty or 6.3-uncertanty. Were the magnitude of the uncertanty is less than one bit of the precision of the number being represented.

To program arround this, the technique used generaly encorporates a value called epsilon. The value of epsilon is chosen such that it is larger than the largest "uncertanty" for the expected numbers to be examined with the IF statement.

Note that the value of uncertanty is different for 6.3 than for 1234.3. If you know what the value ranges are in advance, you can hard wire the epsilon. If not, then prior to the IF you can compute the epsilon, or compute a number approximately equal to your desired epsilon,using the Fortran function EPSILON(x). Where x is a variable of the same type as those to be examined in the IF test.

When you have an epsilon

if((A .gt. (Target - epsilon)) .and. (A .lt. (Target + epsilon)) call FoundA()

Jim Dempsey

[Steven_L_Intel1]

One can also use:

if (abs(A-Target) < SPACING(Target)) ...

SPACING gives you the distance between the argument and the next representable value. You might want to use 2.0*SPACING to be "fuzzier". The Fortran EPSILON intrinsic is not useful for this as it always returns the spacing near 1.0.

[gelarimer]

Thanks JimSteve. Very helpful information. I have been searching on the web for any information/discussions that would help explain aspects of Real Numbers and programming as you have above. Dave Eklund's article 'The Perils of Real Numbers' was a good starting point, and some class room lectures (found on the web)on the IEEE Real Number standard shed some light on that subject, butwould like to geta deeper insight. Is there an in-depth text or reference that would cover the subject from A to Z?

[Steven_L_Intel1]

I haven't run across such myself but I'm not exposed to such things as much as some. I'm sure that if you asked this in comp.lang.fortran you'd get some good pointers. A brief web search turns up a similarly-named article The Perils of Floating Point which I often see cited. When I was in college, I learned a lot from Volume 2 of Donald Knuth's "The Art of Computer Programming" (my favorite of the three volumes that had been published.) Perhaps others here will have additional suggestions.

[gelarimer]

Thanks Steve, appreciate the help and references.