special_union_complement_integer_int64 Function

    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