sets_mod.F90 Source File

#include "macros.h"

module sets_mod
    use constants, only: int32, int64, sp, dp
    use util_mod, only: binary_search_int, stop_all
    use sort_mod, only: sort
    use growing_buffers, only: buffer_int32_1D_t, buffer_int64_1D_t
    implicit none
    private
    public :: subset, set, is_set, is_sorted, special_union_complement, disjoint, &
              operator(.cap.), operator(.complement.), operator(.U.), &
              operator(.in.), operator(.notin.), empty_int

    integer :: empty_int(0) = [integer::]

    !> Check if V is sorted.
    interface is_sorted
        module procedure is_sorted_integer_int64
        module procedure is_sorted_integer_int32
        module procedure is_sorted_real_sp
        module procedure is_sorted_real_dp
    end interface

    !> Check if A and B are disjoint.
    interface disjoint
        module procedure disjoint_integer_int64
        module procedure disjoint_integer_int32
    end interface

    !> Create a set out of A
    interface set
        module procedure set_integer_int64
        module procedure set_integer_int32
    end interface

    !> Check if a given array is a set (ordered and unique elements)
    interface is_set
        module procedure is_set_integer_int64
        module procedure is_set_integer_int32
    end interface

    !> Check if A is a subset of B.
    interface subset
        module procedure subset_integer_int64
        module procedure subset_integer_int32
    end interface

    !> Calculate the union A ∪ B
    interface operator(.U.)
        module procedure union_integer_int64
        module procedure union_integer_int32
    end interface

    !> Calculate the intersection A ∩ B
    interface operator(.cap.)
        module procedure intersect_integer_int64
        module procedure intersect_integer_int32
    end interface

    !> Calculate the complement A / B
    interface operator (.complement.)
        module procedure complement_integer_int64
        module procedure complement_integer_int32
    end interface

    !> Check if element is contained in set.
    interface operator(.in.)
        module procedure test_in_integer_int64
        module procedure test_in_integer_int32
    end interface

    !> Check if element is not contained in set.
    interface operator(.notin.)
        module procedure test_not_in_integer_int64
        module procedure test_not_in_integer_int32
    end interface

    !> Specialiced function with assumptions that speed up performance.
    !> Merge B into A and remove values that are in C.
    !> The result can be written with set notation as A ∪ B / C.
    !> Preconditions (not tested!):
    !>      1. C is a subset of A
    !>      2. A and B are disjoint
    !>      3. B and C are disjoint
    !>      4. A, B, and C are sorted.
    !> The result will be sorted.
    interface special_union_complement
        module procedure special_union_complement_integer_int64
        module procedure special_union_complement_integer_int32
    end interface

contains

    pure function is_sorted_integer_int64 (V, ascending) result(res)
        integer (int64), intent(in) :: V(:)
        logical, intent(in), optional :: ascending
        logical :: ascending_
        logical :: res

        integer :: i

if(present(ascending)) then
    ascending_ = ascending
else
    ascending_ = .true.
endif

        res = .true.
        if (ascending_) then
            do i = 1, size(V) - 1
                if (V(i) > V(i + 1)) then
                    res = .false.
                    return
                end if
            end do
        else
            do i = 1, size(V) - 1
                if (V(i) < V(i + 1)) then
                    res = .false.
                    return
                end if
            end do
        end if
    end function is_sorted_integer_int64
    pure function is_sorted_integer_int32 (V, ascending) result(res)
        integer (int32), intent(in) :: V(:)
        logical, intent(in), optional :: ascending
        logical :: ascending_
        logical :: res

        integer :: i

if(present(ascending)) then
    ascending_ = ascending
else
    ascending_ = .true.
endif

        res = .true.
        if (ascending_) then
            do i = 1, size(V) - 1
                if (V(i) > V(i + 1)) then
                    res = .false.
                    return
                end if
            end do
        else
            do i = 1, size(V) - 1
                if (V(i) < V(i + 1)) then
                    res = .false.
                    return
                end if
            end do
        end if
    end function is_sorted_integer_int32
    pure function is_sorted_real_sp (V, ascending) result(res)
        real (sp), intent(in) :: V(:)
        logical, intent(in), optional :: ascending
        logical :: ascending_
        logical :: res

        integer :: i

if(present(ascending)) then
    ascending_ = ascending
else
    ascending_ = .true.
endif

        res = .true.
        if (ascending_) then
            do i = 1, size(V) - 1
                if (V(i) > V(i + 1)) then
                    res = .false.
                    return
                end if
            end do
        else
            do i = 1, size(V) - 1
                if (V(i) < V(i + 1)) then
                    res = .false.
                    return
                end if
            end do
        end if
    end function is_sorted_real_sp
    pure function is_sorted_real_dp (V, ascending) result(res)
        real (dp), intent(in) :: V(:)
        logical, intent(in), optional :: ascending
        logical :: ascending_
        logical :: res

        integer :: i

if(present(ascending)) then
    ascending_ = ascending
else
    ascending_ = .true.
endif

        res = .true.
        if (ascending_) then
            do i = 1, size(V) - 1
                if (V(i) > V(i + 1)) then
                    res = .false.
                    return
                end if
            end do
        else
            do i = 1, size(V) - 1
                if (V(i) < V(i + 1)) then
                    res = .false.
                    return
                end if
            end do
        end if
    end function is_sorted_real_dp

    pure function set_integer_int64 (V) result(res)
        integer (int64), intent(in) :: V(:)
        integer(int64), allocatable :: res(:)
        integer(int64) :: sorted(size(V)), previous
        type(buffer_int64_1D_t) :: buffer
        integer :: i

        sorted = V
        call sort(sorted)

        call buffer%init()
        if (size(sorted) > 0) then
            call buffer%push_back(sorted(1))
            previous = sorted(1)
        end if
        do i = 2, size(sorted)
            if (previous /= sorted(i)) then
                call buffer%push_back(sorted(i))
                previous = sorted(i)
            end if
        end do
        call buffer%dump_reset(res)
    end function set_integer_int64
    pure function set_integer_int32 (V) result(res)
        integer (int32), intent(in) :: V(:)
        integer(int32), allocatable :: res(:)
        integer(int32) :: sorted(size(V)), previous
        type(buffer_int32_1D_t) :: buffer
        integer :: i

        sorted = V
        call sort(sorted)

        call buffer%init()
        if (size(sorted) > 0) then
            call buffer%push_back(sorted(1))
            previous = sorted(1)
        end if
        do i = 2, size(sorted)
            if (previous /= sorted(i)) then
                call buffer%push_back(sorted(i))
                previous = sorted(i)
            end if
        end do
        call buffer%dump_reset(res)
    end function set_integer_int32

    pure function is_set_integer_int64 (V) result(res)
        integer (int64), intent(in) :: V(:)
        logical :: res
        integer(int64) :: previous
        integer :: i

        res = is_sorted(V)
        if (res) then
            if (size(V) > 0) previous = V(1)
            do i = 2, size(V)
                if (V(i) == previous) then
                    res = .false.
                    return
                else
                    previous = V(i)
                end if
            end do
        end if
    end function is_set_integer_int64
    pure function is_set_integer_int32 (V) result(res)
        integer (int32), intent(in) :: V(:)
        logical :: res
        integer(int32) :: previous
        integer :: i

        res = is_sorted(V)
        if (res) then
            if (size(V) > 0) previous = V(1)
            do i = 2, size(V)
                if (V(i) == previous) then
                    res = .false.
                    return
                else
                    previous = V(i)
                end if
            end do
        end if
    end function is_set_integer_int32

    ! check if A and B are disjoint
    pure function disjoint_integer_int64 (A, B) result(res)
        integer (int64), intent(in) :: A(:), B(:)
        logical :: res
        character(*), parameter :: this_routine = "disjoint_integer_int64"

        integer :: i, j

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(A))) then
            call stop_all (this_routine, "Assert fail: is_sorted(A)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(B))) then
            call stop_all (this_routine, "Assert fail: is_sorted(B)")
        end if
    end block
#endif

        res = .true.
        i = 1; j = 1
        do while (i <= size(A) .and. j <= size(B))
            if (A(i) < B(j)) then
                i = i + 1
            else if (A(i) > B(j)) then
                j = j + 1
            else
                res = .false.
                return
            end if
        end do
    end function disjoint_integer_int64
    pure function disjoint_integer_int32 (A, B) result(res)
        integer (int32), intent(in) :: A(:), B(:)
        logical :: res
        character(*), parameter :: this_routine = "disjoint_integer_int32"

        integer :: i, j

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(A))) then
            call stop_all (this_routine, "Assert fail: is_sorted(A)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(B))) then
            call stop_all (this_routine, "Assert fail: is_sorted(B)")
        end if
    end block
#endif

        res = .true.
        i = 1; j = 1
        do while (i <= size(A) .and. j <= size(B))
            if (A(i) < B(j)) then
                i = i + 1
            else if (A(i) > B(j)) then
                j = j + 1
            else
                res = .false.
                return
            end if
        end do
    end function disjoint_integer_int32

    !> Check if A is a subset of B
    pure function subset_integer_int64 (A, B) result(res)
        integer (int64), intent(in) :: A(:), B(:)
        logical :: res
        character(*), parameter :: this_routine = "subset_integer_int64"

        integer :: i, j

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(A))) then
            call stop_all (this_routine, "Assert fail: is_sorted(A)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(B))) then
            call stop_all (this_routine, "Assert fail: is_sorted(B)")
        end if
    end block
#endif

        if (size(A) == 0) then
            res = .true.
            return
        end if

        res = .false.

        if (size(A) > size(B)) return
        if (A(1) < B(1) .or. A(size(A)) > B(size(B))) return

        i = 1; j = 1
        do while (i <= size(A))
            if (j > size(B)) then
                return
            else if (i == size(A) .and. j == size(B) .and. A(i) /= B(j)) then
                return
            else if (A(i) < B(j)) then
                i = i + 1
            else if (A(i) > B(j)) then
                j = j + 1
            else if (A(i) == B(j)) then
                i = i + 1
                j = j + 1
            end if
        end do
        res = .true.
    end function subset_integer_int64
    pure function subset_integer_int32 (A, B) result(res)
        integer (int32), intent(in) :: A(:), B(:)
        logical :: res
        character(*), parameter :: this_routine = "subset_integer_int32"

        integer :: i, j

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(A))) then
            call stop_all (this_routine, "Assert fail: is_sorted(A)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(B))) then
            call stop_all (this_routine, "Assert fail: is_sorted(B)")
        end if
    end block
#endif

        if (size(A) == 0) then
            res = .true.
            return
        end if

        res = .false.

        if (size(A) > size(B)) return
        if (A(1) < B(1) .or. A(size(A)) > B(size(B))) return

        i = 1; j = 1
        do while (i <= size(A))
            if (j > size(B)) then
                return
            else if (i == size(A) .and. j == size(B) .and. A(i) /= B(j)) then
                return
            else if (A(i) < B(j)) then
                i = i + 1
            else if (A(i) > B(j)) then
                j = j + 1
            else if (A(i) == B(j)) then
                i = i + 1
                j = j + 1
            end if
        end do
        res = .true.
    end function subset_integer_int32

    pure function special_union_complement_integer_int64 (A, B, C) result(D)
        integer (int64), intent(in) :: A(:), B(:), C(:)
        integer (int64) :: D(size(A) + size(B) - size(C))
        character(*), parameter :: this_routine = 'union_complement'

        integer :: i, j, k, l

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(A))) then
            call stop_all (this_routine, "Assert fail: is_sorted(A)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(B))) then
            call stop_all (this_routine, "Assert fail: is_sorted(B)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(C))) then
            call stop_all (this_routine, "Assert fail: is_sorted(C)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (disjoint(A, B))) then
            call stop_all (this_routine, "Assert fail: disjoint(A, B)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (disjoint(B, C))) then
            call stop_all (this_routine, "Assert fail: disjoint(B, C)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (subset(C, A))) then
            call stop_all (this_routine, "Assert fail: subset(C, A)")
        end if
    end block
#endif

        i = 1; j = 1; k = 1; l = 1
        do while (l <= size(D))
            ! Only indices from C have to be added to A
            ! We use assumption that B is a subset of A
            if (i > size(A)) then
                D(l) = B(k)
                k = k + 1
                l = l + 1
                ! Neither indices from B have to be deleted in A
                ! nor indices from C have to be added from C to A.
            else if (j > size(C) .and. k > size(B)) then
                D(l) = A(i)
                i = i + 1
                l = l + 1
                ! No more indices from B have to be deleted in A
            else if (j > size(C)) then
                if (A(i) < B(k)) then
                    D(l) = A(i)
                    i = i + 1
                    l = l + 1
                else if (A(i) > B(k)) then
                    D(l) = B(k)
                    k = k + 1
                    l = l + 1
                end if
                ! No more indices have to be added from C to A
            else if (k > size(B)) then
                if (A(i) /= C(j)) then
                    D(l) = A(i)
                    i = i + 1
                    l = l + 1
                else
                    i = i + 1
                    j = j + 1
                end if
                ! Normal case:
                ! Merge C sorted into A excluding values from B.
            else if (A(i) < B(k)) then
                if (A(i) /= C(j)) then
                    D(l) = A(i)
                    i = i + 1
                    l = l + 1
                else
                    i = i + 1
                    j = j + 1
                end if
            else if (A(i) > B(k)) then
                if (A(i) /= C(j)) then
                    D(l) = B(k)
                    k = k + 1
                    l = l + 1
                else
                    D(l) = B(k)
                    i = i + 1
                    j = j + 1
                    k = k + 1
                    l = l + 1
                end if
            end if
        end do
    end function special_union_complement_integer_int64
    pure function special_union_complement_integer_int32 (A, B, C) result(D)
        integer (int32), intent(in) :: A(:), B(:), C(:)
        integer (int32) :: D(size(A) + size(B) - size(C))
        character(*), parameter :: this_routine = 'union_complement'

        integer :: i, j, k, l

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(A))) then
            call stop_all (this_routine, "Assert fail: is_sorted(A)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(B))) then
            call stop_all (this_routine, "Assert fail: is_sorted(B)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(C))) then
            call stop_all (this_routine, "Assert fail: is_sorted(C)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (disjoint(A, B))) then
            call stop_all (this_routine, "Assert fail: disjoint(A, B)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (disjoint(B, C))) then
            call stop_all (this_routine, "Assert fail: disjoint(B, C)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (subset(C, A))) then
            call stop_all (this_routine, "Assert fail: subset(C, A)")
        end if
    end block
#endif

        i = 1; j = 1; k = 1; l = 1
        do while (l <= size(D))
            ! Only indices from C have to be added to A
            ! We use assumption that B is a subset of A
            if (i > size(A)) then
                D(l) = B(k)
                k = k + 1
                l = l + 1
                ! Neither indices from B have to be deleted in A
                ! nor indices from C have to be added from C to A.
            else if (j > size(C) .and. k > size(B)) then
                D(l) = A(i)
                i = i + 1
                l = l + 1
                ! No more indices from B have to be deleted in A
            else if (j > size(C)) then
                if (A(i) < B(k)) then
                    D(l) = A(i)
                    i = i + 1
                    l = l + 1
                else if (A(i) > B(k)) then
                    D(l) = B(k)
                    k = k + 1
                    l = l + 1
                end if
                ! No more indices have to be added from C to A
            else if (k > size(B)) then
                if (A(i) /= C(j)) then
                    D(l) = A(i)
                    i = i + 1
                    l = l + 1
                else
                    i = i + 1
                    j = j + 1
                end if
                ! Normal case:
                ! Merge C sorted into A excluding values from B.
            else if (A(i) < B(k)) then
                if (A(i) /= C(j)) then
                    D(l) = A(i)
                    i = i + 1
                    l = l + 1
                else
                    i = i + 1
                    j = j + 1
                end if
            else if (A(i) > B(k)) then
                if (A(i) /= C(j)) then
                    D(l) = B(k)
                    k = k + 1
                    l = l + 1
                else
                    D(l) = B(k)
                    i = i + 1
                    j = j + 1
                    k = k + 1
                    l = l + 1
                end if
            end if
        end do
    end function special_union_complement_integer_int32

    !> Return A ∪ B
    !> Assume:
    !>      1. A and B are sorted.
    !> The result will be sorted.
    pure function union_integer_int64 (A, B) result(D)
        integer (int64), intent(in) :: A(:), B(:)
        integer (int64), allocatable :: D(:)
        character(*), parameter :: this_routine = 'union_complement'

        integer (int64), allocatable :: tmp(:)
        integer :: i, j, l

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(A))) then
            call stop_all (this_routine, "Assert fail: is_sorted(A)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(B))) then
            call stop_all (this_routine, "Assert fail: is_sorted(B)")
        end if
    end block
#endif

        allocate(tmp(size(A) + size(B)))

        i = 1; j = 1; l = 1
        do while (l <= size(tmp))
            if (i > size(A) .and. j > size(B)) then
                exit
            else if (i > size(A)) then
                tmp(l) = B(j)
                j = j + 1
                l = l + 1
            else if (j > size(B)) then
                tmp(l) = A(i)
                i = i + 1
                l = l + 1
            else if (A(i) == B(j)) then
                tmp(l) = A(i)
                i = i + 1
                j = j + 1
                l = l + 1
            else if (A(i) < B(j)) then
                tmp(l) = A(i)
                i = i + 1
                l = l + 1
            else if (A(i) > B(j)) then
                tmp(l) = B(j)
                j = j + 1
                l = l + 1
            end if
        end do

        associate(final_size => l - 1)
            D = tmp(:final_size)
        end associate
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(D))) then
            call stop_all (this_routine, "Assert fail: is_sorted(D)")
        end if
    end block
#endif
    end function union_integer_int64
    pure function union_integer_int32 (A, B) result(D)
        integer (int32), intent(in) :: A(:), B(:)
        integer (int32), allocatable :: D(:)
        character(*), parameter :: this_routine = 'union_complement'

        integer (int32), allocatable :: tmp(:)
        integer :: i, j, l

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(A))) then
            call stop_all (this_routine, "Assert fail: is_sorted(A)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(B))) then
            call stop_all (this_routine, "Assert fail: is_sorted(B)")
        end if
    end block
#endif

        allocate(tmp(size(A) + size(B)))

        i = 1; j = 1; l = 1
        do while (l <= size(tmp))
            if (i > size(A) .and. j > size(B)) then
                exit
            else if (i > size(A)) then
                tmp(l) = B(j)
                j = j + 1
                l = l + 1
            else if (j > size(B)) then
                tmp(l) = A(i)
                i = i + 1
                l = l + 1
            else if (A(i) == B(j)) then
                tmp(l) = A(i)
                i = i + 1
                j = j + 1
                l = l + 1
            else if (A(i) < B(j)) then
                tmp(l) = A(i)
                i = i + 1
                l = l + 1
            else if (A(i) > B(j)) then
                tmp(l) = B(j)
                j = j + 1
                l = l + 1
            end if
        end do

        associate(final_size => l - 1)
            D = tmp(:final_size)
        end associate
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(D))) then
            call stop_all (this_routine, "Assert fail: is_sorted(D)")
        end if
    end block
#endif
    end function union_integer_int32

    !> Return A ∩ B
    !> Assume:
    !>      1. A and B are sorted.
    !> The result will be sorted.
    pure function intersect_integer_int64 (A, B) result(D)
        integer (int64), intent(in) :: A(:), B(:)
        integer (int64), allocatable :: D(:)
        character(*), parameter :: this_routine = 'union_complement'

        integer (int64), allocatable :: tmp(:)
        integer :: i, j, l

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(A))) then
            call stop_all (this_routine, "Assert fail: is_sorted(A)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(B))) then
            call stop_all (this_routine, "Assert fail: is_sorted(B)")
        end if
    end block
#endif

        allocate(tmp(min(size(A), size(B))))

        i = 1; j = 1; l = 1
        do while (l <= size(tmp))
            if (i > size(A) .or. j > size(B)) then
                exit
            else if (A(i) == B(j)) then
                tmp(l) = A(i)
                i = i + 1
                j = j + 1
                l = l + 1
            else if (A(i) < B(j)) then
                i = i + 1
            else if (A(i) > B(j)) then
                j = j + 1
            end if
        end do

        associate(final_size => l - 1)
            D = tmp(:final_size)
        end associate
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(D))) then
            call stop_all (this_routine, "Assert fail: is_sorted(D)")
        end if
    end block
#endif
    end function intersect_integer_int64
    pure function intersect_integer_int32 (A, B) result(D)
        integer (int32), intent(in) :: A(:), B(:)
        integer (int32), allocatable :: D(:)
        character(*), parameter :: this_routine = 'union_complement'

        integer (int32), allocatable :: tmp(:)
        integer :: i, j, l

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(A))) then
            call stop_all (this_routine, "Assert fail: is_sorted(A)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(B))) then
            call stop_all (this_routine, "Assert fail: is_sorted(B)")
        end if
    end block
#endif

        allocate(tmp(min(size(A), size(B))))

        i = 1; j = 1; l = 1
        do while (l <= size(tmp))
            if (i > size(A) .or. j > size(B)) then
                exit
            else if (A(i) == B(j)) then
                tmp(l) = A(i)
                i = i + 1
                j = j + 1
                l = l + 1
            else if (A(i) < B(j)) then
                i = i + 1
            else if (A(i) > B(j)) then
                j = j + 1
            end if
        end do

        associate(final_size => l - 1)
            D = tmp(:final_size)
        end associate
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(D))) then
            call stop_all (this_routine, "Assert fail: is_sorted(D)")
        end if
    end block
#endif
    end function intersect_integer_int32

    !> Return A / B
    !> Assume:
    !>      1. A and B are sorted.
    !> The result will be sorted.
    pure function complement_integer_int64 (A, B) result(D)
        integer (int64), intent(in) :: A(:), B(:)
        integer (int64), allocatable :: D(:)
        character(*), parameter :: this_routine = 'union_complement'

        integer (int64), allocatable :: tmp(:)
        integer :: i, j, l

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(A))) then
            call stop_all (this_routine, "Assert fail: is_sorted(A)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(B))) then
            call stop_all (this_routine, "Assert fail: is_sorted(B)")
        end if
    end block
#endif

        allocate(tmp(size(A)))

        i = 1; j = 1; l = 1
        do while (l <= size(tmp))
            if (i > size(A)) then
                exit
            else if (j > size(B)) then
                tmp(l) = A(i)
                i = i + 1
                l = l + 1
            else if (A(i) == B(j)) then
                i = i + 1
                j = j + 1
            else if (A(i) < B(j)) then
                tmp(l) = A(i)
                i = i + 1
                l = l + 1
            else if (A(i) > B(j)) then
                j = j + 1
            end if
        end do

        associate(final_size => l - 1)
            D = tmp(:final_size)
        end associate
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(D))) then
            call stop_all (this_routine, "Assert fail: is_sorted(D)")
        end if
    end block
#endif
    end function complement_integer_int64
    pure function complement_integer_int32 (A, B) result(D)
        integer (int32), intent(in) :: A(:), B(:)
        integer (int32), allocatable :: D(:)
        character(*), parameter :: this_routine = 'union_complement'

        integer (int32), allocatable :: tmp(:)
        integer :: i, j, l

#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(A))) then
            call stop_all (this_routine, "Assert fail: is_sorted(A)")
        end if
    end block
#endif
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(B))) then
            call stop_all (this_routine, "Assert fail: is_sorted(B)")
        end if
    end block
#endif

        allocate(tmp(size(A)))

        i = 1; j = 1; l = 1
        do while (l <= size(tmp))
            if (i > size(A)) then
                exit
            else if (j > size(B)) then
                tmp(l) = A(i)
                i = i + 1
                l = l + 1
            else if (A(i) == B(j)) then
                i = i + 1
                j = j + 1
            else if (A(i) < B(j)) then
                tmp(l) = A(i)
                i = i + 1
                l = l + 1
            else if (A(i) > B(j)) then
                j = j + 1
            end if
        end do

        associate(final_size => l - 1)
            D = tmp(:final_size)
        end associate
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(D))) then
            call stop_all (this_routine, "Assert fail: is_sorted(D)")
        end if
    end block
#endif
    end function complement_integer_int32

    pure function test_in_integer_int64 (element, set) result(res)
        integer (int64), intent(in) :: element, set(:)
        character(*), parameter :: this_routine = 'test_in_integer_int64'
        logical :: res
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(set))) then
            call stop_all (this_routine, "Assert fail: is_sorted(set)")
        end if
    end block
#endif
        res = binary_search_int(set, element) /= -1
    end function
    pure function test_in_integer_int32 (element, set) result(res)
        integer (int32), intent(in) :: element, set(:)
        character(*), parameter :: this_routine = 'test_in_integer_int32'
        logical :: res
#ifdef DEBUG_
    block
        use util_mod, only: stop_all
        if (.not. (is_sorted(set))) then
            call stop_all (this_routine, "Assert fail: is_sorted(set)")
        end if
    end block
#endif
        res = binary_search_int(set, element) /= -1
    end function

    pure function test_not_in_integer_int64 (element, set) result(res)
        integer (int64), intent(in) :: element, set(:)
        logical :: res
        res = .not. (element .in. set)
    end function
    pure function test_not_in_integer_int32 (element, set) result(res)
        integer (int32), intent(in) :: element, set(:)
        logical :: res
        res = .not. (element .in. set)
    end function

end module sets_mod