splitmix64_next_poisson_int Function

   function splitmix64_next_poisson_int(gen, mean_milli, err) result(r)
      !! Poisson count with mean `mean_milli/1000`, from a splitmix64_t generator.
      type(splitmix64_t), intent(inout) :: gen
         !! Generator, advanced once per inter-arrival drawn.
      integer(default_int), intent(in) :: mean_milli
         !! A thousand times the mean, in `[0, POISSON_MAX_MEAN_MILLI]`.
      type(error_t), intent(inout), optional :: err
         !! Set to `ERROR_VALIDATION` when `mean_milli` is out of range.
      integer(int32) :: r

      integer(int64) :: gap_mean, window, elapsed

      r = 0_int32
      if (mean_milli < 0_default_int .or. mean_milli > POISSON_MAX_MEAN_MILLI) then
         if (present(err)) then
            call err%set(ERROR_VALIDATION, &
                         "pic_random_dist: next_poisson_int mean_milli out of range")
         end if
         return
      end if
      if (mean_milli == 0_default_int) return

      ! Mean gap between arrivals, in window ticks: window / (mean_milli/1000).
      ! That division truncates, so the gap actually used is a little short of
      ! the one asked for, and the count over a fixed window comes out high --
      ! by up to 0.048%, worst at mean_milli = 999596.
      !
      ! Fixing the gap and deriving the window from it instead turns that
      ! round-off to the sampler's advantage: whatever gap the truncation
      ! landed on, the window is the one for which that gap is exactly right.
      ! The residual is the truncation of the window itself, under 1e-6
      ! relative.
      gap_mean = (POISSON_WINDOW*1000_int64)/int(mean_milli, int64)
      window = (gap_mean*int(mean_milli, int64))/1000_int64
      elapsed = 0_int64
      do
         elapsed = elapsed + splitmix64_next_exponential_int(gen, gap_mean)
         if (elapsed > window) exit
         r = r + 1_int32
      end do
   end function splitmix64_next_poisson_int