YARP - Yet another rounding problem

[h-faber]

Hi all,

consider the following lines of code:

Code:

REAL RUNDST,FAKT
         RUNDST = 0.3000
         DO 52, I=1,10
              IF (I/(RUNDST*10) .EQ. 1.) THEN
                   FAKT=0.1**I
                   FAKT=(FAKT /0.1**I) /10**I
                   GOTO 54
              END IF
52       CONTINUE

You might guess it - the if-clause is never fulfilled so the loop never runs into the THEN-part. Although the condition should fit for I=3. VS Debugger anyway shows I/(RUNDST*10) as 0.9999999.

I have no idea why 3/(0.3*10) should not equal 1.0 as the division should be performed *after* the 0.3 turns into 3.0. Even using some REAL variable TMP=I does not help.

Is there any known solution for this problem? Any help appreciated.

[Steven_L_Intel1]

0.1 (and 0.3) are not exactly representable in binary floating point, so the values you actually get are slightly higher or lower, depending on the value. 0.1 ends up slightly lower, so adding or multiplying it just increases the difference.

[h-faber]

---

sblionel wrote:
0.1 (and 0.3) are not exactly representable in binary floating point, so the values you actually get are slightly higher or lower, depending on the value. 0.1 ends up slightly lower, so adding or multiplying it just increases the difference.

---

Hi Steve,

thanks for the explanation, but I would like to have a solution for this problem. :smileywink:
I think this is a common problem and I am not the only one looking for a solution. Anyway I also haven't found some intrinsic ROUND function. Which advice can you give, how can I get an "equal" for my condition?

[TimP]

NINT() and ANINT() have been part of Fortran for 27 years. You're entitled to make your own ROUND() function if the name is important.

[Steven_L_Intel1]

What is usually done in such cases is compare the absolute value of the difference against some small value deemed as "small enough".

[h-faber]

---

tim18 wrote:
NINT() and ANINT() have been part of Fortran for 27 years. You're entitled to make your own ROUND() function if the name is important.

---

NINT() does not work for my value of approx. 156.3625 and also turns it into 156.362 (when val*1000, NINT(val), val/1000)
Maybe ANINT does the job, will try it out later.

[h-faber]

---

sblionel wrote:
What is usually done in such cases is compare the absolute value of the difference against some small value deemed as "small enough".

---

Yes this is probably the way I will go.

[Steven_L_Intel1]

Neither NINT nor ANINT will work for you. Both round from .5 towards the nearest integer, but as you note, that's not what you want. ANINT has a real datatype, so you get 3.0 and not 3, which NINT gives.

[jim_dempsey]

Many decimal fractions require an infinite series of binary bits to represent. You seeasimilareffect on a decimal calculator when performing 1/3 (0.3333333333) which is not 1/3. When comparing values that were derrived from fractions you must take this in to consideration.

One way is to use an epsilon value.

epsilon = TINY(value)

...

if(value .ge. target-epsilon .and. value .le. target+epsilon) call foo

The other method is to recode such that fractions are eliminated in the computations that are destined for use in an IFstatement. Example

Change:

TargetSeconds = 0.3

...

if(ValueInSeconds .eq. TargetSeconds) call foo

to

TargetMilliseconds = 300.

...

if(ValueInMilliseconds .eq. TargetMilliseconds) call foo

[h-faber]

Hi Jim,

my task was a bit more flexible as here is a function which does the rounding and takes some arguments, also the position where to round. 0.3 for thousandth, 0.2 for hundredth etc. So for me it is precise enough to divide this 0.3 or so by 1000. and check if the rounded value differs more or less than my 0.001 x rounding precision.

Code:

         DO 52, I=1,10
              TEMP = I/(RUNDST*10)
              IF (ABS(TEMP - 1.) .LE. (RUNDST*DELTAFAKT)) THEN
C              DO SOMETHING
                   GOTO 54
              END IF
52       CONTINUE
        ENDIF
54    CONTINUE

where RUNDST is 0.2, 0.3 etc. like described above and DELTAFAKT is 0.001 by hard.

[ketan9]

Use the following code, it works for me

REAL RUNDST,FAKT
RUNDST = 3.0D-01
DO 52, I=1,10
IF (I/INT(RUNDST*1.0D01) .EQ. 1) THEN
FAKT=0.1**I
FAKT=(FAKT /0.1**I) /10**I
GOTO 54
END IF
52 CONTINUE

Note: You must try declaring real variables using D... eg.
Real tmp
Instead of declaring like
tmp=0.001
declare it as
tmp=1.0D-03
By using this declaration syntax, you would not have any trailing numbers.
In previous case, tmp would be declared as tmp=0.00100000009 something like this where as in later case, it would be declared exactly as you want 0.001

Hope this solves your problem:)
Ketan

[Intel_C_Intel]

Hello.

Compiling with /real-size:64 /fpconstant /QxN and defining all constants with full double precision is required to gain maximum reliability of floating point calculations:

DOUBLE PRECISION PI
PARAMETER (PI = 3.14159265358979323846D00)

Best Regards,

Lars Petter Endresen