      module RandomMethods

      use iso_fortran_env

      !the following subroutine is part of the DCDFLIB library acquired
      !from http://www.netlib.org/random/
      !description: various TOMS algorithms and Abramowitz&Stegun;
      !             Bux Dekker zero-finding
      !authors: Barry W. Brown, James Lovato,
      !         Kathy Russell <bwb@odin.mda.uth.tmc.edu>
      !for compiling issues of both the library and the program, see:
      !http://stanford.edu/class/me200c/tutorial_77/18.3_libraries.html
      external cdfchi !Chi-Square function; see DCDFLIB documentation

      integer(kind=int32) :: XorX32 = 314159265;
      logical :: showHistograms_ = .true.

      contains

      !the following subroutine is taken from the following website:
      !http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html
      subroutine init_random_seed()
        implicit none
        integer, allocatable :: seed(:)
        integer :: i, n, un, istat, dt(8), pid, t(2), s
        integer(8) :: count, tms

        call random_seed(size = n)
        allocate(seed(n))
        ! First try if the OS provides a random number generator
        open(newunit=un, file="/dev/urandom", access="stream",
     &       form="unformatted", action="read", status="old",
     &       iostat=istat)
        if (istat == 0) then
           read(un) seed
           close(un)
        else
           ! Fallback to XOR:ing the current time and pid. The PID is
           ! useful in case one launches multiple instances of the same
           ! program in parallel.
           call system_clock(count)
           if (count /= 0) then
              t = transfer(count, t)
           else
              call date_and_time(values=dt)
              tms = (dt(1) - 1970) * 365_8 * 24 * 60 * 60 * 1000
     &            + dt(2) * 31_8 * 24 * 60 * 60 * 1000
     &            + dt(3) * 24 * 60 * 60 * 60 * 1000
     &            + dt(5) * 60 * 60 * 1000
     &            + dt(6) * 60 * 1000 + dt(7) * 1000
     &            + dt(8)
              t = transfer(tms, t)
           end if
           s = ieor(t(1), t(2))
           pid = getpid() + 1099279 ! Add a prime
           s = ieor(s, pid)
           if (n >= 3) then
              seed(1) = t(1) + 36269
              seed(2) = t(2) + 72551
              seed(3) = pid
              if (n > 3) then
                 seed(4:) = s + 37 * (/ (i, i = 0, n - 4) /)
              end if
           else
              seed = s + 37 * (/ (i, i = 0, n - 1 ) /)
           end if
        end if
        call random_seed(put=seed)
      end subroutine init_random_seed

      integer*8 function LCG(modulus,multiplier,increment,seed,modified)
      !Linear congruential generator
      !  0 <  modulus
      !  0 <  multiplier < modulus
      !  0 <= increment  < modulus
      !  0 <= seed       < modulus
      !The modified parameter is a custom modification, turned off by
      !default. If set true, time in milliseconds will be used to alter
      !seed.
        integer*8 modulus,multiplier,increment,seed
        logical, optional :: modified
        logical alter

        integer date_time(8)
        integer*8 add_seed

        add_seed=int(0,8)
        alter=.false.
        if (present(modified)) alter=modified
        if (alter) then
          call date_and_time(values=date_time)
          if (alter .and. increment+date_time(8)<modulus)
     &      add_seed=increment+date_time(8)
        end if

        LCG=MOD(multiplier*(seed+add_seed)+increment,modulus)
        return
      end function

      subroutine runLCG(modified)
      !run Linear Congruential Generator
      !The modified parameter is a custom modification, turned off by
      !default. If set true, time in milliseconds will be used to alter
      !seed.
        integer*8 LCGm,LCGa,LCGc,LCGx0,LCGtotal,LCGi,LCGx
        integer order
        logical continue,alter
        logical, optional :: modified
        character(1) choice
        character(5) OutFormat
        integer*8, dimension(:), allocatable :: arrHisto

100     format(A)

        alter=.false.
        if (present(modified)) alter=modified

        print*
        print*
        print*,"Linear Congruential Generator (LCG)"
        print*,"-----------------------------------"
        print*
        if (alter) print*,"NOTE that time (in ms) will be used for "//
     &                    "seed modification."
        if (alter) print*
        print*,"This generator is defined by the equation:"
        print*
        print*,"    X[n+1] := (a*X[n] + c) mod m"
        print*
        print*,"where"
        print*,'  m is the "modulus":                  0 < m'
        print*,'  a is the "multiplier":               0 < a < m'
        print*,'  c is the "increment":                0 <= c < m'
        print*,'  X[0] is the start value ("seed"):    0 <= X[0] < m'
        print*
        print*,"Please enter these values to initiate the random "//
     &         "number generator!"
        call InputI8(LCGm, "  m = ",int(0,8))
        call InputI8(LCGa, "  a = ",int(0,8),LCGm)
        call InputI8(LCGc, "  c = ",int(-1,8),LCGm)
        call InputI8(LCGx0,"  X[0] = ",int(-1,8),LCGm)
        print*
        if (allocated(arrHisto)) deallocate(arrHisto)
        allocate(arrHisto(LCGm))
        arrHisto(:)=int(0,8)

        call InputI8(LCGtotal,"How many numbers would you like to "//
     &                          "have generated? ",int(0,8))
        print*

        order=0
        do i=0,11
          order=i
          if (LCGm/10**i<1) exit
        end do
        write(OutFormat,"(A,I2,A)")"(I",order+2,")"

        LCGi=int(0,8)
        LCGx=LCGx0
        continue=alter
        if (alter) print*,"Calculation started..."
        do
          !23 lines per circle of computation
          do i=1,23
             if (.not. continue) write(*,"(A)",advance="no")" "
            !caculate, write into histogram, output as many numbers as
            !possible each line
            do j=1,78/(order+2)
              LCGi=LCGi+int(1,8)
              !do not continue if total amount is reached
              if (LCGi>LCGtotal) exit
              LCGx=LCG(LCGm,LCGa,LCGc,LCGx,alter)
              arrHisto(LCGx+1)=arrHisto(LCGx+1)+int(1,8)
              if (.not. continue) write(*,OutFormat,advance="no")LCGx
            end do
            if (.not. continue) print*
            if (LCGi>LCGtotal) exit
          end do
          if (.not. continue) print*
          if (LCGi>LCGtotal) exit
          if (.not. continue) then
            choice=""
210         select case (choice)
              case ("d")
                print*
              case ("c")
                print*
                print*,"Calculation underway..."
                continue=.true.
              case default
                write(*,100,advance="no")" Would you like to "//
     &                                   "[d]isplay the next page "//
     &                                   "or [c]ontinue silently? "
                read*,choice
                goto 210
            end select
          end if
        end do
        if (continue) print*

        print*,"Calculation complete."
        print*

        call Analyze("Linear Congruential Generator",arrHisto,
     &               showHistograms_)

        deallocate(arrHisto)
        return
      end subroutine

      subroutine runXorShift32()
      !run XorShift32() function
        integer*8 XorSmax,currXorS
        integer*8 XorSseed,XorStotal,XorSi
        integer order
        logical continue
        character(1) choice
        character(5) OutFormat
        integer*8, dimension(:), allocatable :: arrHisto

100     format(A)

        print*
        print*
        print*,"Xorshift Random Number Generator"
        print*,"--------------------------------"
        print*
        print*,"Please enter upper boundary of random number "//
     &         "distribution to be created!"
        call InputI8(XorSmax,"  max = ",int(0,8))
        print*
        print*,"Please enter a non-zero seed value to initiate the "//
     &         "random number generator!"
        call InputI8(XorSseed,"  seed = ",int(0,8))

        print*
        if (allocated(arrHisto)) deallocate(arrHisto)
        allocate(arrHisto(XorSmax+1))
        arrHisto(:)=int(0,8)

        call InputI8(XorStotal,"How many numbers would you like to "//
     &                         "have generated? ",int(0,8))
        print*

        order=0
        do i=0,11
          order=i
          if (XorSmax/10**i<1) exit
        end do
        write(OutFormat,"(A,I2,A)")"(I",order+2,")"

        XorSi=int(0,8)
        continue=.false.
        call init_XorShift32(int(XorSseed,int32))
        do
          !23 lines per circle of computation
          do i=1,23
             if (.not. continue) write(*,"(A)",advance="no")" "
            !caculate, write into histogram, output as many numbers as
            !possible each line
            do j=1,78/(order+2)
              XorSi=XorSi+int(1,8)
              !do not continue if total amount is reached
              if (XorSi>XorStotal) exit
              currXorS=anint(XorShift32()*(XorSmax+1)-0.5d0,8)
              arrHisto(currXorS+1)=arrHisto(currXorS+1)+int(1,8)
              if (.not. continue)
     &          write(*,OutFormat,advance="no")currXorS
            end do
            if (.not. continue) print*
            if (XorSi>XorStotal) exit
          end do
          if (.not. continue) print*
          if (XorSi>XorStotal) exit
          if (.not. continue) then
            choice=""
210         select case (choice)
              case ("d")
                print*
              case ("c")
                print*
                print*,"Calculation underway..."
                continue=.true.
              case default
                write(*,100,advance="no")" Would you like to "//
     &                                   "[d]isplay the next page "//
     &                                   "or [c]ontinue silently? "
                read*,choice
                goto 210
            end select
          end if
        end do
        if (continue) print*

        print*,"Calculation complete."
        print*

        call Analyze("Xorshift Random Number Generator",arrHisto,
     &               showHistograms_)

        deallocate(arrHisto)
        return
      end subroutine

      subroutine runRAND()
      !run Fortran's built-in RAND function
        integer*8 RANDmax,currRAND,RANDtotal,RANDi
        integer RANDseed,order
        logical continue
        character(1) choice
        character(5) OutFormat
        integer*8, dimension(:), allocatable :: arrHisto

100     format(A)

        print*
        print*
        print*,"Fortran's built-in RAND function"
        print*,"--------------------------------"
        print*
        print*,"Please enter upper boundary of random number "//
     &         "distribution to be created!"
        call InputI8(RANDmax,"  max = ",int(0,8))
        print*
        print*,"Please enter a seed value to initiate the random "//
     &         "number generator!"
        call InputI(RANDseed,"  seed = ",-1)

        print*
        if (allocated(arrHisto)) deallocate(arrHisto)
        allocate(arrHisto(RANDmax+1))
        arrHisto(:)=int(0,8)

        call InputI8(RANDtotal,"How many numbers would you like to "//
     &                         "have generated? ",int(0,8))
        print*

        order=0
        do i=0,11
          order=i
          if (RANDmax/10**i<1) exit
        end do
        write(OutFormat,"(A,I2,A)")"(I",order+2,")"

        RANDi=int(0,8)
        continue=.false.
        call srand(RANDseed)
        do
          !23 lines per circle of computation
          do i=1,23
             if (.not. continue) write(*,"(A)",advance="no")" "
            !caculate, write into histogram, output as many numbers as
            !possible each line
            do j=1,78/(order+2)
              RANDi=RANDi+int(1,8)
              !do not continue if total amount is reached
              if (RANDi>RANDtotal) exit
              currRAND=anint(rand(0)*(RANDmax+1)-0.5,8)
              arrHisto(currRAND+1)=arrHisto(currRAND+1)+int(1,8)
              if (.not. continue)
     &          write(*,OutFormat,advance="no")currRAND
            end do
            if (.not. continue) print*
            if (RANDi>RANDtotal) exit
          end do
          if (.not. continue) print*
          if (RANDi>RANDtotal) exit
          if (.not. continue) then
            choice=""
210         select case (choice)
              case ("d")
                print*
              case ("c")
                print*
                print*,"Calculation underway..."
                continue=.true.
              case default
                write(*,100,advance="no")" Would you like to "//
     &                                   "[d]isplay the next page "//
     &                                   "or [c]ontinue silently? "
                read*,choice
                goto 210
            end select
          end if
        end do
        if (continue) print*

        print*,"Calculation complete."
        print*

        call Analyze("Fortran''s built-in RAND function",arrHisto,
     &               showHistograms_)

        deallocate(arrHisto)
        return
      end subroutine

      subroutine runRANDOM_NUMBER()
      !run Fortran's built-in RANDOM_NUMBER function
        double precision RANDreal
        integer*8 RANDseed,RANDmax,currRAND,RANDtotal,RANDi
        integer order
        logical continue
        character(1) choice
        character(5) OutFormat
        integer*8, dimension(:), allocatable :: arrHisto

100     format(A)

        print*
        print*
        print*,"Fortran's built-in RANDOM_NUMBER function"
        print*,"-----------------------------------------"
        print*
        print*,"NOTE that the program's process ID and time (in ms) "//
     &         "will be used as seed."
        print*
        print*,"Please enter upper boundary of random number "//
     &         "distribution to be created!"
        call InputI8(RANDmax,"  max = ",int(0,8))

        print*
        if (allocated(arrHisto)) deallocate(arrHisto)
        allocate(arrHisto(RANDmax+1))
        arrHisto(:)=int(0,8)

        call InputI8(RANDtotal,"How many numbers would you like to "//
     &                         "have generated? ",int(0,8))
        print*

        order=0
        do i=0,11
          order=i
          if (RANDmax/10**i<1) exit
        end do
        write(OutFormat,"(A,I2,A)")"(I",order+2,")"

        RANDi=int(0,8)
        continue=.false.
        call init_random_seed()
        do
          !23 lines per circle of computation
          do i=1,23
             if (.not. continue) write(*,"(A)",advance="no")" "
            !caculate, write into histogram, output as many numbers as
            !possible each line
            do j=1,78/(order+2)
              RANDi=RANDi+int(1,8)
              !do not continue if total amount is reached
              if (RANDi>RANDtotal) exit
              call random_number(RANDreal)
              currRAND=anint(RANDreal*(RANDmax+1)-0.5,8)
              arrHisto(currRAND+1)=arrHisto(currRAND+1)+int(1,8)
              if (.not. continue)
     &          write(*,OutFormat,advance="no")currRAND
            end do
            if (.not. continue) print*
            if (RANDi>RANDtotal) exit
          end do
          if (.not. continue) print*
          if (RANDi>RANDtotal) exit
          if (.not. continue) then
            choice=""
210         select case (choice)
              case ("d")
                print*
              case ("c")
                print*
                print*,"Calculation underway..."
                continue=.true.
              case default
                write(*,100,advance="no")" Would you like to "//
     &                                   "[d]isplay the next page "//
     &                                   "or [c]ontinue silently? "
                read*,choice
                goto 210
            end select
          end if
        end do
        if (continue) print*

        print*,"Calculation complete."
        print*

        call Analyze("Fortran''s built-in RANDOM_NUMBER function",
     &               arrHisto,showHistograms_)

        deallocate(arrHisto)
        return
      end subroutine

      subroutine Analyze(title,arrHisto,showHistograms)
      !analizes random numbers in arrHisto; shows histogram and count
      !distribution if so defined in showHistograms (which is an
      !optional parameter and .true. by default)
        character(*) title
        integer*8, dimension(:), allocatable :: arrHisto,arrDistri
        logical, optional :: showHistograms
        logical showHist
        integer*8 maxDistri,totalRNs
        double precision chiSq,chiSqNexp,chiSqP,chiSqQ,chiSqBound,
     &                   chiSqDF,KSS0,KSdmax5,KSdumax,KSdomax,KSF0,KSdu,
     &                   KSdo,KScurrN
        integer chiSqStatus
        character(6) tResult

        showHist=.true.
        if (present(showHistograms)) showHist=showHistograms

        if (showHist) then
          !show histogram
          print*,"Showing histogram..."
          call showHistogram("RNG_Histo","Random Number Histogram",
     &                       title,arrHisto)
          print*

          !distribution: what is the maximal count any number has?
          maxDistri=int(0,8)
          do i=1,size(arrHisto)
            if (arrHisto(i)>maxDistri) maxDistri=arrHisto(i)
          end do

          !allocate distribution array and fill with values
          allocate(arrDistri(maxDistri+1))
          arrDistri(:)=int(0,8)
          do i=1,size(arrHisto)
            arrDistri(arrHisto(i)+1)=arrDistri(arrHisto(i)+1)+int(1,8)
          end do

          !show histogram distribution
          print*,"Showing histogram distribution..."
          call showHistogram("RNG_Distr",
     &                       "Random Number Histogram Distribution",
     &                       title,arrDistri,.true.)
          deallocate(arrDistri)
          print*
        end if !showHist

        print*
        print*,"(The following tests assume equally distributed "//
     &         "random numbers.)"
        print*

        !general calculations
        totalRNs=0
        do i=1,size(arrHisto)
          totalRNs=totalRNs+arrHisto(i)
        end do

        !Chi squared test
        tResult="failed"
        chiSq=0.d0
        chiSqNexp=1.d0*totalRNs/size(arrHisto)
        do i=1,size(arrHisto)
          chiSq=chiSq+(arrHisto(i)-chiSqNexp)**2/chiSqNexp
        end do
        chiSqDF=size(arrHisto)*1.d0
        call cdfchi(1,chiSqP,chiSqQ,chiSq,chiSqDF,chiSqStatus,
     &              chiSqBound)
        if (chiSqQ<0.05) tResult="passed"

        print*,"Chi-square test:"
        write(*,"(A25,F20.4)")"computed Chi squared:",chiSq
        write(*,"(A25,F20.4)")"level of significance:",chiSqQ
        write(*,"(A25,A20)")"result for alpha=0.05:",tResult
        print*

        !Kolmogorow-Smirnow test
        if (totalRNs<40) then
          print*,"Kolmogorow-Smirnow test not performed (N too low)."
          print*
        else
          tResult="passed"
          KSS0=0.d0
          !KSdmax5=1.36d0/dsqrt(totalRNs*1.d0)
          !d_max for alpha=0.05
          KSdmax5=1.d0*dsqrt(log(2.d0/0.05d0))/dsqrt(2.d0*totalRNs)
          KSdumax=0.d0
          KSdomax=0.d0
          KScurrN=0.d0
          do i=1,size(arrHisto)
            KSF0=0.d0
            KScurrN=KScurrN+1.d0*arrHisto(i)
            KSF0=1.d0*KScurrN/totalRNs
            KSdu=dabs(KSS0-KSF0) !for d_ui, S(x_i-1) is used
            KSS0=1.d0*i/size(arrHisto)
            KSdo=dabs(KSS0-KSF0) !for d_oi, S(x) is used
            if (KSdu>KSdumax) KSdumax=KSdu
            if (KSdo>KSdomax) KSdomax=KSdo
          end do
          if (KSdumax>KSdmax5) tResult="failed"
          if (KSdomax>KSdmax5) tResult="failed"

          print*,"Kolmogorow-Smirnow test:"
          write(*,"(A25,F20.4)")"greatest abs. distance:",
     &                          dmax1(KSdumax,KSdomax)
          write(*,"(A25,F20.4)")"should be lower than:",KSdmax5
          write(*,"(A25,A20)")"result for alpha=0.05:",tresult
          print*
        end if

        return
      end subroutine

      subroutine ShowHistogram(filename,title,description,arrHisto,
     &                         autozoom)
      !saves arrHisto to file and plots it using gnuplot
      !parameter autozoom neglects empty columns (zeros) left and right
      !of the occurance of data in arrHisto; it is .false. by default
      !HINT: it is taken into account that arrHisto(1) refers to the
      !      value at 0 as arrHisto(0) is not defined - that is to say:
      !      this particular shift is taken care of
        character(*) filename,title,description
        logical, optional :: autozoom
        logical zoom
        integer*8, dimension(:), allocatable :: arrHisto
        integer*8 xrange,yrange,xstart

100     format(A)
101     format(I10,I10)
102     format(A,I12,A)
103     format(A,I12,A,I12,A)

        zoom=.false.
        if (present(autozoom)) zoom=autozoom

        xrange=size(arrHisto)-int(1,8)
        xstart=int(0,8)
        if (zoom) then
          do i=1,size(arrHisto)
            if (arrHisto(i)>0) exit
            xstart=i-int(1,8)
          end do
          do i=size(arrHisto),1,-1
            if (arrHisto(i)>0) exit
            xrange=i-int(1,8)
          end do
        end if

        yrange=int(1,8)
        do i=1,size(arrHisto)
          if (arrHisto(i)>yrange) yrange=arrHisto(i)
        end do
        yrange=ceiling(yrange*1.2)

        open(1,file=filename//".dat",position='rewind')
        write(1,100)"#"//title
        write(1,100)'#Histogram data: '//description
        !when accessing the histogram array, note that 0 is a possible
        !result, but is saved in arrHisto(1)
        do i=1,xrange+1
          write(1,101),i-1,arrHisto(i)
        end do
        close(1)

        open(1,file=filename//".plt",position='rewind')
        write(1,100)"set title '"//title//"';"
        write(1,100)"set grid;"
        write(1,103)"set xrange[",xstart,":",xrange,"];"
        write(1,102)"set yrange[0:",yrange,"];"
        write(1,100)"set xlabel '';"
        write(1,100)"set ylabel 'N';"
        !write(1,100)"set logscale y;"
        !write(1,100)"set format y '10^%L';"
        write(1,100)"plot '"//filename//".dat' title' "//description//
     &              "' with boxes;"
        write(1,100) "pause mouse key"
        write(1,100) "if (MOUSE_KEY != 27) reread"
        close(1)

        call system("gnuplot "//filename//".plt")
        return
      end subroutine

      subroutine InputI(Variable,Text,minValue,maxValue)
      !Input procedure with check for boundaries; Variable will be set
      !HINT: boundaries are NOT included, i.e. minValue<Varible, not <=
      !      and maxValue>Variable, not >=
        integer Variable,lower,upper
        character(*) Text
        integer, optional :: minValue,maxValue

        !use lowest and highest possible numbers as default boundaries
        lower=1-huge(0_2)
        if (present(minValue)) lower=minValue
        upper=huge(0_2)
        if (present(maxValue)) upper=maxValue
        if (lower>upper) then
          print*,"ERROR: Lower boundary greater than upper boundary."
          stop
        end if

        do
          write(*,"(A)",advance="no")" "//Text
          read*,Variable
          if (lower<Variable .and. Variable<upper) exit
          print*,"Value entered is incorrect. Respect boundaries!"
        end do
        return
      end subroutine

      subroutine InputI8(Variable,Text,minValue,maxValue)
      !Input procedure with check for boundaries; Variable will be set
      !HINT: boundaries are NOT included, i.e. minValue<Varible, not <=
      !      and maxValue>Variable, not >=
        integer*8 Variable,lower,upper
        character(*) Text
        integer*8, optional :: minValue,maxValue

        !use lowest and highest possible numbers as default boundaries
        lower=1-huge(0_4)
        if (present(minValue)) lower=minValue
        upper=huge(0_4)
        if (present(maxValue)) upper=maxValue
        if (lower>upper) then
          print*,"ERROR: Lower boundary greater than upper boundary."
          stop
        end if

        do
          write(*,"(A)",advance="no")" "//Text
          read*,Variable
          if (lower<Variable .and. Variable<upper) exit
          print*,"Value entered is incorrect. Respect boundaries!"
        end do
        return
      end subroutine

      subroutine init_XorShift32(seed)
      !initializes the XorShift32 random number generator; if the
      !optional parameter seed is not set or zero, 314159265 is used
      integer(kind=int32), optional :: seed
        XorX32=314159265;
        if (present(seed) .and. seed/=0) XorX32=seed;
      end subroutine

      double precision function XorShift32()
      !actually, XorShift uses unsigned integers; as Fortran does not
      !allow this, a double precision value in [0,1] is assigned
        !do actual XorShift operation
        XorX32 = xor(XorX32,ishft(XorX32,+13))
        XorX32 = xor(XorX32,ishft(XorX32,-17))
        XorX32 = xor(XorX32,ishft(XorX32,+5))

        !transform to double precision
        if (XorX32>0) XorShift32=1.d0*XorX32/huge(i8)
        if (XorX32<0) XorShift32=XorX32/(-1.d0*huge(i8)-1.d0)
        return
      end function

      end module


      program Random

      use RandomMethods

      character(1) choice

100   format(A)

      print*,"RANDOM NUMBER GENERATORS"
      print*,"========================"

      choice=""
200   select case (choice)
        case ("a")
          !Linear Congruential Generator
          call runLCG(.false.)
        case ("b")
          !Linear Congruential Generator, modified
          call runLCG(.true.)
        case ("c")
          !Xorshift Random Number Generator
          call runXorShift32()
        case ("d")
          !Fortran's built-in RAND function
          call runRAND()
        case ("e")
          !Fortran's built-in RANDOM_NUMBER function
          call runRANDOM_NUMBER()
        case ("x")
          !exit
          goto 999
        case default
          print*
          print*,"Please choose the method you'd like to analyze:"
          print*
          print*,"  [a] Linear Congruential Generator (LCG)"
          print*,"  [b] LCG using time (in ms) as seed modification"
          print*,"  [c] Xorshift Random Number Generator"
          print*,"  [d] Fortran's built-in RAND function"
          print*,"  [e] Fortran's built-in RANDOM_NUMBER function"
          print*
          print*,"  [x] Exit"
          print*
          write(*,100,advance="no")" Your choice: "
          read*,choice
          goto 200
      end select

      print*,"Returning to main menu..."
      print*
      choice=""
      goto 200

999   print*
      print*,"Done."

      end
