CODEWITHSUNDEEP
fortran
mattiav27
mattiav27

Reputation: 685

Small number treated as zero when adding in Fortran

I am writing a code for a Monte Carlo simulation in Fortran, but I am having a lot of problems because of the small numbers involved.

The biggest problem is that in my code particle positions are not updated; the incriminated code looks like this

x=x+step*cos(p)*sin(t)

with step=0.001. With this, the code won't update the position and I get a infinite loop because the particle never exits the region. If I modify my code with something like this:

 x=x+step

or 

x=x+step*cos(t)

there is no problem. So it seems that the product step*cos(t)*cos(p)(of the order 10**-4) is too small and is treated as zero.

x is of the order 10**4.

How do I solve this problem in portable way?

My compiler is the latest f95.

Upvotes: 1

Views: 773

Answers (1)

francescalus
francescalus

Reputation: 32366

Your problem is essentially the one of this other question. However, it's useful to add some Fortran-specific comments.

As in that other question, the discrete nature of floating point numbers mean that there is a point where one number is too small to make a difference when added to another. In the case of this question:

if (1e4+1e-4==1e4) print *, "Oh?"
if (1d4+1d-4==1d4) print *, "Really?"
end

That is, you may be able to use double precision reals and you'll see the problem go away.

What is the smallest number you can add to 1e4 to get something different from 1e4 (or to 1d4)?

print *, 1e4 + SPACING(1e4), 1e4+SPACING(1e4)/2
print *, 1d4 + SPACING(1d4), 1d4+SPACING(1d4)/2
end

This spacing varies with the size of the number. For large numbers it is large and around 1 it is small.

print*, EPSILON(1e0), SPACING([(1e2**i,i=0,5)])
print*, EPSILON(1d0), SPACING([(1d2**i,i=0,5)])
end

Upvotes: 2

Related Questions