Search code examples
randomscopefortranfortran77

Linear congruential generator - output is all 0s?


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.


Solution

  • 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!).
    • The subroutine 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.
    • I used 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...
    • I named the program file as "test.f90" (note the suffix .f90), which allows free-format. Also, it is okay to use lower-case letters for variables.
    • [Also, because 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 :)