! **************************************************** !
! ******************* lcg-cpp11.f90 ****************** !
! **************************************************** !
!
!    Leonardo Dagdug, Ivan Pompa-García, Jason Peña
!
! From the book:
! **************************************************** !
!      Diffusion Under Confinement:
!               A Journey Through Counterintuition    
! **************************************************** !

! **************************************************** !
! This file contains the source code that implements
! the minstd function specified in the cpp11
! standard. It generates 10,000 numbers and compares
! the result with the one stated in the standard.
! **************************************************** !

! **************************************************** !
! To compile with GFortran, you must include the
! helpers.f90 module.
!         gfortran helpers.f90 lcg_cpp11.f90
! Then you can run the program:
!         ./a.out
!
! After running, a file called «minstd_numbers.dat»
! is generated, and can be plotted with any software
! of your choice. To do this using gnuplot, execute:
!         gnuplot -p -e 'plot "minstd_numbers.dat"'
! -p means that the plot must persist until the
!    user closes it explicitly.
! -e is used to specify the instructions
!    for using  gnuplot directly in the commandline.
! **************************************************** !

program lcg_cpp11

! Load the helpers module, which contains functions,
! constants, and more...
use helpers

! Every data type must be explicitly stated.
implicit none

! Declare and initialize our LCG parameters
! as constants.
integer(kind=i32), parameter :: a = 0
integer(kind=i32), parameter :: b = 48271
! We must write the result. If we try to calculate it,
! the compiler drops an error because 2^31 overflows the
! capacity of an integer, even if 2^31 - 1 fits.
integer(kind=i32), parameter :: m = 2147483647

! How many numbers will be generated.
integer(kind=i32), parameter :: n = 10000
! According to c++11 standard, the 10,000th generated 
! number must be 399268537 when the initial seed is 1.
integer(kind=i64), parameter :: xstd = 399268537

! A unit number for the file where the numbers
! will be saved. Fortran >= 2003 will
! assign this automatically in the open call.
integer(kind=i32) :: FUNIT

! Number in the sequence. Initially this must be
! the seed.
integer(kind=i64) :: x = 1

! A counter
integer :: i

! Open a file to write the normalized numbers
open(newunit=FUNIT, file='minstd_numbers.dat')

do i=1, n
  x = a + mod(b * x, m)

  ! You must not forget that «x», and «m» are integers
  ! and Fortran will be perform integer arithmetic.
  ! You must multiply one of the operands by a real
  ! factor.
  write(FUNIT,*) 1.0_dp * x / m
end do

! Close the file
close(unit=FUNIT)

if(x == xstd) then
  print *, 'The 10000 number from seed 1 is correct: ' &
    // nstr(x)
else
  print *, &
    'Something goes wrong, the last number should be ' &
    // nstr(xstd) // ' and not ' // nstr(x)
end if

end program lcg_cpp11