constrained_sample_fast Subroutine

    subroutine constrained_sample_fast(this, contain, contain_ilut, renormalization, pos, val, prob)
        class(AliasSampler_t), intent(in) :: this
        integer, intent(in) :: contain(:)
        integer(n_int), intent(in) :: contain_ilut(0 : )
        real(dp), intent(in) :: renormalization
        integer, intent(out) :: pos, val
        real(dp), intent(out) :: prob
        routine_name("constrained_sample_fast")
        ASSERT(is_set(contain))
        ASSERT(1 <= contain(1) .and. contain(size(contain)) <= size(this%probs%ptr))
        ASSERT(renormalization .isclose. (sum(this%get_prob(contain))))
        ASSERT(size(contain) == sum(popcnt(contain_ilut)))

        if (near_zero(renormalization)) then
            pos = 0; val = 0; prob = 1.0
            return
        end if

        if (renormalization < redrawing_cutoff) then
        block
            type(CDF_Sampler_t) :: cdf_sampler
            cdf_sampler = CDF_Sampler_t(this%probs%ptr(contain), renormalization)
            call cdf_sampler%sample(pos, prob)
            val = contain(pos)
            return
        end block
        end if

        val = this%table%sample()
        do while (.not. IsOcc(contain_ilut, val))
            val = this%table%sample()
        end do
        pos = int(binary_search_int(contain, val))
        prob = this%probs%ptr(val) / renormalization
        ASSERT(prob .isclose. (this%probs%ptr(val) / sum(this%probs%ptr(contain))))
    end subroutine constrained_sample_fast