mtprng_rand64 Function

    function mtprng_rand64(state) result(r)

        !! Obtain the next 32-bit integer in the psuedo-random sequence
        !! Uses the Mersenne Twister algorithm

        type(mtprng_state), intent(inout) :: state
        integer(INT64) :: r

        ! internal constants
        integer(INT64), dimension(0:1), parameter :: mag01 = [ 0_INT64, -1727483681_INT64 ]

        ! Period parameters
        integer(INT64), parameter :: UPPER_MASK =  2147483648_INT64
        integer(INT64), parameter :: LOWER_MASK =  2147483647_INT64

        ! Tempering parameters
        integer(INT64), parameter :: TEMPERING_B = -1658038656_INT64
        integer(INT64), parameter :: TEMPERING_C =  -272236544_INT64

        ! Note: variable names match those in original example
        integer(INT32) :: kk

        ! Generate N words at a time
        if (state%mti >= mtprng_N) then
            ! The value -1 acts as a flag saying that the seed has not been set.
            if (state%mti == -1) call mtprng_init(4357_INT32,state)

            ! Fill the mt array
            do kk = 0, mtprng_N - mtprng_M - 1
                r = ior(iand(state%mt(kk),UPPER_MASK),iand(state%mt(kk+1),LOWER_MASK))
                state%mt(kk) = ieor(ieor(state%mt(kk + mtprng_M),ishft(r,-1_INT64)),mag01(iand(r,1_INT64)))
            end do

            do kk = mtprng_N - mtprng_M, mtprng_N - 2
                r = ior(iand(state%mt(kk),UPPER_MASK),iand(state%mt(kk+1),LOWER_MASK))
                state%mt(kk) = ieor(ieor(state%mt(kk + (mtprng_M - mtprng_N)),ishft(r,-1_INT64)),mag01(iand(r,1_INT64)))
            end do

            r = ior(iand(state%mt(mtprng_N-1),UPPER_MASK),iand(state%mt(0),LOWER_MASK))
            state%mt(mtprng_N-1) = ieor(ieor(state%mt(mtprng_M-1),ishft(r,-1)),mag01(iand(r,1_INT64)))

            ! Start using the array from first element
            state%mti = 0
        end if

        ! Here is where we actually calculate the number with a series of
        !   transformations
        r = state%mt(state%mti)
        state%mti = state%mti + 1

        r = ieor(r,ishft(r,-11))
        r = iand(4294967295_INT64,ieor(r,iand(ishft(r, 7),TEMPERING_B)))
        r = iand(4294967295_INT64,ieor(r,iand(ishft(r,15),TEMPERING_C)))
        r = ieor(r,ishft(r,-18))

    end function mtprng_rand64