forums.silverfrost.com Forum Index forums.silverfrost.com
Welcome to the Silverfrost forums
 
 FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

output problem

 
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Plato
View previous topic :: View next topic  
Author Message
dzovan137



Joined: 12 Sep 2013
Posts: 4

PostPosted: Thu Sep 12, 2013 4:08 pm    Post subject: output problem Reply with quote

I have been having trouble with following code:

Code:

PROGRAM HitMiss
integer :: hits,i,N
real :: r,s,p

! RANDOM STUFF FROM THE CLOCK JUST TO MAKE SURE THAT WITH EVERY RUN THE INITIAL CONDITIONS ARE A LITTLE BIT DIFFERENT
INTEGER ::  gj, clock
INTEGER, DIMENSION(:), ALLOCATABLE :: seed
CALL RANDOM_SEED(size = gj)
ALLOCATE(seed(gj))       
CALL SYSTEM_CLOCK(COUNT=clock)
seed = clock + 37 * (/ (i - 1, i = 1, gj) /)
CALL RANDOM_SEED(PUT = seed)
DEALLOCATE(seed)

!Trial number
N = 50000

hits = 0
do i=1,N
 
 call random_number(r)
 call random_number(s)

r = r +1
s = s*2**10 -1


if ( s <= r**10 -1) then
  hits = hits + 1

  end if

 

end do

! Here there is an error in the compiler because I get zero which is obviously wrong
p =(hits/N)*(2**10 -1)
write (*,*) 'Hits:', hits
write (*,*) 'Integral:',p



END PROGRAM HitMiss



The final output with:

Code:

write (*,*) 'Integral:',p


gives me zero which shouldn't be the case.

Can someone help me with this?

Cheerz,

Jovan
Back to top
View user's profile Send private message
PaulLaidler
Site Admin


Joined: 21 Feb 2005
Posts: 7916
Location: Salford, UK

PostPosted: Fri Sep 13, 2013 6:44 am    Post subject: Reply with quote

Try using the debugger to step through and see where the calculation goes wrong.
Back to top
View user's profile Send private message AIM Address
LitusSaxonicum



Joined: 23 Aug 2005
Posts: 2388
Location: Yateley, Hants, UK

PostPosted: Fri Sep 13, 2013 11:25 am    Post subject: Reply with quote

There are 3 possible cases comparing hits to N:

hits < N

hits = N

and hits > N (but given your DO loop, this is impossible)

Without running your problem, it seems that the most likely outcome is the first, i.e. hits < N, in which case, in integer arithmetic, hits/N = 0 , and flowing directly from that, p=0.0.

If you force real arithmetic, then p won't be zero.
For example:

Code:
TOP = hits
BOTTOM = N
P = (TOP/BOTTOM)*(etc)


The compiler compiled your code exactly as the you wrote it.

Eddie
Back to top
View user's profile Send private message
JohnCampbell



Joined: 16 Feb 2006
Posts: 2554
Location: Sydney

PostPosted: Fri Sep 13, 2013 1:32 pm    Post subject: Reply with quote

Try using:
p =(real(hits)/real(N))*(2**10 -1)

John
Back to top
View user's profile Send private message
LitusSaxonicum



Joined: 23 Aug 2005
Posts: 2388
Location: Yateley, Hants, UK

PostPosted: Fri Sep 13, 2013 3:08 pm    Post subject: Reply with quote

Hi John,

REAL converts to the default real precision, whereas an assignment converts to the precision of the target, I think. This saves me needing to think about REAL, REAL with a KIND specifier, FLOAT, DBLE and all the rest (usually extensions)!

Eddie
Back to top
View user's profile Send private message
brucebowler
Guest





PostPosted: Fri Sep 13, 2013 3:57 pm    Post subject: Reply with quote

But... When does the TYPE promotion occur (my feeble mind can't recall)? If it occurs AFTER the computation, then INT/INT will result in a truncated INT, not a REAL (of what ever precision). Of course if the TYPE promotion occurs before the assignment then John's REALs will have no effect (but do no harm, either).
Back to top
LitusSaxonicum



Joined: 23 Aug 2005
Posts: 2388
Location: Yateley, Hants, UK

PostPosted: Fri Sep 13, 2013 5:10 pm    Post subject: Reply with quote

If the INT/INT is in brackets, then the subexpression is computed in integer mode. I use this all the time, and I am confident that is the case. The 2^10-1 is certainly done in integer, because that is bracketed, and in any case, is all integer. The type promotion to real (p) occurs with the assignment.

As done by John, the type promotions for hits and N are done in the function calls, with the last bracketed subexpression done in integer. Then, with a real*integer ... I'd expect that all to be done in real. This is convenient, because the assignment is to real. If p happened to be integer, I would then expect real computation followed by a truncation.

The point of the assignments to temporary variables is that one then never needs to wonder where the type promotion occurs, or which function to use (I can never remember, and don't expect anyone else to).

There's a chance that you only have to convert hits OR N, and not both, to get the result wanted, but doing both is belt and braces.

Eddie
Back to top
View user's profile Send private message
dzovan137



Joined: 12 Sep 2013
Posts: 4

PostPosted: Sun Sep 15, 2013 4:31 pm    Post subject: Re: Reply with quote

LitusSaxonicum wrote:
There are 3 possible cases comparing hits to N:

hits < N

hits = N

and hits > N (but given your DO loop, this is impossible)

Without running your problem, it seems that the most likely outcome is the first, i.e. hits < N, in which case, in integer arithmetic, hits/N = 0 , and flowing directly from that, p=0.0.

If you force real arithmetic, then p won't be zero.
For example:

Code:
TOP = hits
BOTTOM = N
P = (TOP/BOTTOM)*(etc)


The compiler compiled your code exactly as the you wrote it.

Eddie



The reassigment of values to a temp variables didn't work.

But redefining
Code:
real :: TOP,BOTTOM
made it worked.
Back to top
View user's profile Send private message
dzovan137



Joined: 12 Sep 2013
Posts: 4

PostPosted: Sun Sep 15, 2013 4:37 pm    Post subject: Re: Reply with quote

JohnCampbell wrote:
Try using:
p =(real(hits)/real(N))*(2**10 -1)

John


This does the job. I think there is no loss by forcing the type conversion in the case of my code.
Back to top
View user's profile Send private message
LitusSaxonicum



Joined: 23 Aug 2005
Posts: 2388
Location: Yateley, Hants, UK

PostPosted: Sun Sep 15, 2013 4:52 pm    Post subject: Reply with quote

Well it worked for me with the 3 lines I gave you.

As your code is presented, you can simply insert TOP and BOTTOM because they are implicitly defined as REAL type, and while modern programming styles insist on the definition of types for all variables, the default in Fortran is as it always was.

If the code you supplied was a cut down version for brevity, then the original could have contained extra lines, such as an IMPLICIT NONE statement. If it had, then TOP and BOTTOM would have been of unknown type, and there would have been a compilation error. Any IMPLICIT statement or type statement that caused TOP or BOTTOM to be INTEGER would have stopped it working. Many people who post on this forum would advocate the use of IMPLICIT NONE.


The initial problem remains that you did an INTEGER/INTEGER divide, with the numerator less than the denominator, and the result of that is zero.

Eddie
Back to top
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> Plato All times are GMT + 1 Hour
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © 2001, 2005 phpBB Group