Fortran - Having to add integer command

  • Fortran
  • Thread starter Nugso
  • Start date
  • #1
Nugso
Gold Member
170
10

Main Question or Discussion Point

Hello everyone. Today I learnt how to write a program which calculates and prints integer squares roots from 1 to 1000. Here is what I wrote:
Code:
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;

Code:
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.
 

Answers and Replies

  • #2
jtbell
Mentor
15,578
3,557
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.
 
  • #3
11,872
5,526
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
 
  • #5
11,872
5,526
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
 
  • #6
Nugso
Gold Member
170
10
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.
 

Related Threads on Fortran - Having to add integer command

  • Last Post
Replies
4
Views
4K
  • Last Post
Replies
7
Views
5K
  • Last Post
Replies
7
Views
1K
  • Last Post
Replies
4
Views
2K
Replies
3
Views
2K
Replies
6
Views
14K
Replies
3
Views
1K
Replies
7
Views
4K
Replies
1
Views
3K
  • Last Post
Replies
2
Views
2K
Top