I've been trying to make a pretty basic LCG pseudorandom number generator in Fortran 77 to print out 1000 random numbers to a file, but for whatever reason the output is just 1000 0s. The entire code is pretty short so I've combed it over multiple times and tried changing some things around but I can't for the life of me figure out what is wrong. I have a hunch that it could be a scope issue (if such a concept is even useful in Fortran), but that is really unfounded.
PROGRAM RANDOM
COMMON ISEED, RANDOMNUMBER
ISEED = 123
OPEN (UNIT=1,FILE='rand.in',STATUS='UNKNOWN')
J=1
7 CALL RANDU(ISEED)
J=J+1
WRITE(1,*) RANDOMNUMBER
IF(J<1000)GOTO 7
STOP
END
SUBROUTINE RANDU(ISEED)
PARAMETER (IMAX = 2147483647, IMAXINV = 1./IMAX)
ISEED = ISEED * 65539
IF(ISEED<0) ISEED = ISEED + IMAX + 1
RANDOMNUMBER = ISEED * IMAXINV
RETURN
END
Does anyone have any ideas here? I'm fresh out.
OK now to augment the answer by @Jerry101, I have written a modified code. Here, the key problem is that IMAXINV
is not explicitly declared as REAL
, so it is interpreted as INTEGER
(as a result, IMAXINV = 1.0 / IMAX
becomes always 0 in the original code). Also, I have removed ISEED
from the COMMON
block (because it is passed as an argument) and put another COMMON
statement in RANDU
to share variables among routines. With these modifications the program seems to work correctly.
PROGRAM RANDOM
COMMON RANDOMNUMBER !<--- ISEED is deleted from here
ISEED = 123
J=1
7 CALL RANDU(ISEED)
J=J+1
WRITE(*,*) RANDOMNUMBER !<--- write to STDOUT for test
IF (J < 100) GOTO 7
END
SUBROUTINE RANDU(ISEED)
real IMAXINV !<--- this is necessary
COMMON RANDOMNUMBER !<--- this is also necessary to share variables
PARAMETER (IMAX = 2147483647, IMAXINV = 1./IMAX)
ISEED = ISEED * 65539
IF(ISEED<0) ISEED = ISEED + IMAX + 1
RANDOMNUMBER = ISEED * IMAXINV
END
As suggested in the other Answer, we could also use FUNCTION
to return a variable directly. Then we do not need to use COMMON
, so the code becomes a bit cleaner.
PROGRAM RANDOM
ISEED = 123
J=1
7 RANDOMNUMBER = RANDU(ISEED)
J=J+1
WRITE(*,*) RANDOMNUMBER
IF (J < 100) GOTO 7
END
FUNCTION RANDU(ISEED)
real IMAXINV
PARAMETER (IMAX = 2147483647, IMAXINV = 1./IMAX)
ISEED = ISEED * 65539
IF(ISEED<0) ISEED = ISEED + IMAX + 1
RANDU = ISEED * IMAXINV !<--- "RANDU" is the return variable
END
But note that when FUNCTION
is used, the type of return variable should be declared explicitly in the calling routine if the function name does not conform to the implicit rule. (In the above code, RANDU
is not declared explicitly because it is interpreted as REAL
). So anyway, there are many caveats in the implicit typing rule in Fortran77...
Additional notes:
To avoid these pitfalls, I suggest using Fortran >=90 (rather than Fortran77) because it provides many features for preventing such errors. For example, a minimally modified code may look like this:
module mymodule
contains
subroutine randu ( istate, ran )
implicit none
integer, parameter :: IMAX = 2147483647
real, parameter :: IMAXINV = 1.0 / IMAX
integer, intent(inout) :: istate
real, intent(out) :: ran
istate = istate * 65539
if ( istate < 0 ) istate = istate + IMAX + 1
ran = istate * IMAXINV
end subroutine
end module
program main
use mymodule, only: randu
implicit none
integer :: j, istate
real :: randomnumber
istate = 123 !! seed for RANDU()
do j = 1, 99
call randu ( istate, randomnumber )
write(*,*) randomnumber
enddo
end program
Here,
implicit none
is used to enforce the declaration of all variables explicitly. This is useful to help avoid incorrect typing of variables (such as IMAXINV
in the Question!).RANDU
is contained in a module
so that the compiler provides explicit interface and many useful checking (in short, module
is something like a namespace in C++). module
can also be used to define global variables in a way much safer than COMMON
.do
... enddo
construct for looping over j
rather than incrementing it manually and using goto
. The former is actually easier to use, and also goto
tends to make a code often less readable...iseed
stores the information on the current status of (pseudo) random number generator, it may be better to use some different varaible name (like istate etc?) to remind that its value needs to be kept during calls.]So if you are interested, please consider using a more modern version of Fortran (rather than Fortran77) which allows us to write safer and more robust codes :)