Fortran - Having to add integer command

[Nugso]

Hello everyone. Today I learned how to write a program which calculates and prints integer squares roots from 1 to 1000. Here is what I wrote:

program example
real::i
do i=1,1000
if(sqrt(i)**2==i) then
print*,i
end if
end do
end program example

After compling this, 25, 49 etcs can't be seen. It lists like;
1.0000
4.0000
9.0000

and so on. However, if I write it like this;

program example
real::i
integer::physicsforum
do i=1,1000
physicsforum=sqrt(i)
if(physicsforum**2=i)then
print*,i
end if
end do
end program example

then I do see 25.0000, 49.0000 etc too. What is confusing me is why do I have to add integer for sqrt(i) to see 25.0000, 49.0000?

Thanks.

 

[jtbell]

You should never compare 'real' numbers for exact equality, because arithmetic inevitably introduces small errors because of the finite number of binary bits used. Your first program probably compares 25.00000 with 25.00001 or 24.99999 or something like that.

When you need to do a comparision like this, you must allow for small deviations from equality by doing something like this:

if (abs(a - b) < epsilon)

where 'epsilon' is a suitably small number.

Or you can try to set up the comparision so that it uses integers instead, like you did in the second program. Here you need to be careful because when you assign a real value to an integer variable, it is truncated to the nearest lower integer, not rounded (which could go either up or down). So 24.99999 would become 24, not 25.

 

[jedishrfu]

By assigning it to an integer you've thrown away the small decimal remainder that is common to floating pt numbers.

This is why programmers never do equals testing like you did but instead do this

if ( abs(sqrt(i)**2 - i) < 0.000001) ...

to determine if they are equal to within some small range

 

[Nugso]

Thank you very much for your real quick replies. But why are only 25 and 49 are not shown? I'm uploading the picture of compiled-program within a minute.

http://i.imgur.com/kvKw5Rk.png
http://i.imgur.com/fEgY2Vk.png

 

[jedishrfu]

because they have some stray bits that the other values don't have.

you'd need to print the values out with maximal decimal formatting to see what the stray bits may be.

For example you might print 1.000,000 but the number is really 1.000,000,001 (commas added for clarity) in memory.

There are also even smaller guard bits that aren't shown normally but can be teased out be multiplying the number up a few orders of magnitude up.

see wikipedia for something on guard bits:

http://en.wikipedia.org/wiki/Guard_digit

 

[Nugso]

because they have some stray bits that the other values don't have.

you'd need to print the values out with maximal decimal formatting to see what the stray bits may be.

For example you might print 1.000,000 but the number is really 1.000,000,001 (commas added for clarity) in memory.

There are also even smaller guard bits that aren't shown normally but can be teased out be multiplying the number up a few orders of magnitude up.

see wikipedia for something on guard bits:

http://en.wikipedia.org/wiki/Guard_digit

Thank you very much sir.