! gfortran -Ofast -C -march=native -mtune=native -mavx512f -g -Wall -Wextra -fimplicit-none -malign-double -funroll-loops -ftree-vectorize -fopt-info -ftree-vectorizer-verbose=1 -fcheck=all -fsanitize=address threefry_pi_estimate.f90 -o threefry_pi_estimate && time perf stat ./threefry_pi_estimate
include 'threefry.f90'

program threefry_pi_estimate
    use iso_fortran_env, only: int64, pr => real64
    use ThreefryModule, only: threefry, threefry_key, threefry_e, threefry_result
    implicit none

    integer (int64), parameter :: SEED = int (z'1234', int64)
    integer (int64), parameter :: total_count = 1600000000
    integer (int64), parameter :: histo_size = 100
    integer (int64), dimension (histo_size, histo_size) :: histo = 0
    real (pr), parameter :: one = 1, pi = acos (-one)
    real (pr), parameter :: MULT = 1.0_pr / real (huge (1_int64), pr)
   !real (pr), dimension (total_count) :: data
    integer (pr), dimension (0:1) :: c
    real (pr), dimension (0:1) :: c_r
    type (threefry_result) :: c_e
    type (threefry_key) :: K
    integer (int64) :: i
    integer :: pi_count = 0
    real (pr) :: pi_estimate

    K%k(0) = 0_int64
    K%k(1) = SEED

    do concurrent (i = 1:total_count)
        ! c = threefry ([i, 0_int64], K)
        ! c_r = c * MULT
        c_e = threefry_e (threefry_result([i, 0_int64]), K)
        c_r = c_e%v * MULT
!        data (i) = c_r (0)
        histo (int (c_r (0) * histo_size + 1), int (c_r (1) * histo_size + 1)) = &
        histo (int (c_r (0) * histo_size + 1), int (c_r (1) * histo_size + 1)) + 1

       !if (DOT_PRODUCT (c_r, c_r) < 1.0_pr) then
        if (hypot (c_r (0), c_r (1)) < 1.0_pr) then
                pi_count = pi_count + 1
        end if
    end do

    pi_estimate = 4 * real (pi_count, pr) / real (total_count, pr)
    print *, "Estimate of pi:", &
             pi_estimate, &
             1 / sqrt (real (total_count, pr)), &
             (pi_estimate / pi - 1)

    print *, "histo:",   &
    sum (histo),    ";", &
    sum (histo, 1) * histo_size / real (total_count, pr), ";", &
    sum (histo, 2) * histo_size / real (total_count, pr), ";", &
    sqrt (sum (histo**2, 1) / real (total_count, pr) - (sum (histo, 1) / real (total_count, pr))**2), ";", &
    sqrt (sum (histo**2, 2) / real (total_count, pr) - (sum (histo, 2) / real (total_count, pr))**2) ! , ";",
end program threefry_pi_estimate