Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Fortran - Having to add integer command

  1. Mar 6, 2013 #1

    Nugso

    User Avatar
    Gold Member

    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 (Text):
    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 (Text):
    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.
     
  2. jcsd
  3. Mar 6, 2013 #2

    jtbell

    User Avatar

    Staff: Mentor

    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.
     
  4. Mar 6, 2013 #3

    jedishrfu

    Staff: Mentor

    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. Mar 6, 2013 #4

    Nugso

    User Avatar
    Gold Member

    Last edited: Mar 6, 2013
  6. Mar 6, 2013 #5

    jedishrfu

    Staff: Mentor

    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
     
  7. Mar 7, 2013 #6

    Nugso

    User Avatar
    Gold Member

    Thank you very much sir.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Fortran - Having to add integer command
  1. Fortran DO command (Replies: 4)

Loading...