From a8f185342ad2461a08e773454775fd48bdeab3ee Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Sat, 9 Sep 2023 10:49:37 +0100
Subject: [PATCH 01/17] support union of multiple intervals

---
 src/interval.jl  | 36 +++++++++++++++++++++++++++++-------
 test/runtests.jl |  4 ++++
 2 files changed, 33 insertions(+), 7 deletions(-)

diff --git a/src/interval.jl b/src/interval.jl
index 569d417..906b241 100644
--- a/src/interval.jl
+++ b/src/interval.jl
@@ -119,14 +119,36 @@ end
 
 intersect(d1::AbstractInterval, d2::AbstractInterval) = intersect(Interval(d1), Interval(d2))
 
+union(I::TypedEndpointsInterval...) = unioninterval(intervalunion!([I...]))
+unioninterval(x) = isone(length(x)) ? x[1] : throw(ArgumentError("Cannot construct union of disjoint sets."))
+function intervalunion!(I::Vector{<:TypedEndpointsInterval})
+    sort!(I, by = leftendpoint)
+
+    # filter out empty sets
+    k = 1
+    while k <= length(I)
+        if isempty(I[k])
+            popat!(I,k)
+        else
+            k += 1
+        end
+    end
 
-function union(d1::TypedEndpointsInterval{L1,R1,T1}, d2::TypedEndpointsInterval{L2,R2,T2}) where {L1,R1,T1,L2,R2,T2}
-    T = promote_type(T1,T2)
-    isempty(d1) && return Interval{L2,R2,T}(d2)
-    isempty(d2) && return Interval{L1,R1,T}(d1)
-    any(∈(d1), endpoints(d2)) || any(∈(d2), endpoints(d1)) ||
-        throw(ArgumentError("Cannot construct union of disjoint sets."))
-    _union(d1, d2)
+    # merge intervals
+    k = 2
+    while k <= length(I)
+        x = leftendpoint(I[k])
+        y = rightendpoint(I[k-1])
+        if x > y
+            k += 1
+        elseif x < y || x ∈ I[k] || x ∈ I[k-1]
+            I[k-1] = _union(I[k-1],I[k])
+            popat!(I,k)
+        else
+            k += 1
+        end
+    end
+    I
 end
 
 # these assume overlap
diff --git a/test/runtests.jl b/test/runtests.jl
index e36a964..e6be894 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -408,6 +408,10 @@ struct IncompleteInterval <: AbstractInterval{Int} end
             i5 = nextfloat(one(T)) .. 2one(T)
             i_empty = one(T) ..zero(T)
 
+            # - union of multiple intervals
+            # issue #103
+            @test i1 ∪ i2 ∪ i3 ∪ i4 ∪ i5 ∪ i_empty == 0..3
+
             # - union of completely overlapping intervals
             # i1      0 ------>------ 1
             # i2         1/3 ->- 1/2

From f34a57e380ca3037bed830719aaf1fe5470f2a3e Mon Sep 17 00:00:00 2001
From: Tianyi Pu <44583944+putianyi889@users.noreply.github.com>
Date: Tue, 12 Sep 2023 16:21:22 +0100
Subject: [PATCH 02/17] Update test/runtests.jl

Co-authored-by: Yuto Horikawa <hyrodium@gmail.com>
---
 test/runtests.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index e6be894..684688c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -410,7 +410,10 @@ struct IncompleteInterval <: AbstractInterval{Int} end
 
             # - union of multiple intervals
             # issue #103
-            @test i1 ∪ i2 ∪ i3 ∪ i4 ∪ i5 ∪ i_empty == 0..3
+            @test i1 ∪ i3 ∪ i4 == 0..3  # Equivalent to (i1 ∪ i3) ∪ i4
+            @test_throws ArgumentError i1 ∪ i4 ∪ i3  # The order of union matters
+            @test ∪(i1, i3, i4) == 0..3
+            @test ∪(i1, i4, i3) == 0..3
 
             # - union of completely overlapping intervals
             # i1      0 ------>------ 1

From 88cc31cfc109b38d8fc68460bc9055fc0d1e70f3 Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 19 Feb 2024 19:41:59 +0000
Subject: [PATCH 03/17] final setup

---
 Project.toml           |  5 ++-
 src/interval.jl        | 42 ++++++++-----------------
 src/unionalgorithms.jl | 70 ++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 86 insertions(+), 31 deletions(-)
 create mode 100644 src/unionalgorithms.jl

diff --git a/Project.toml b/Project.toml
index 1c7af6f..882f04a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,10 +1,12 @@
 name = "IntervalSets"
 uuid = "8197267c-284f-5f27-9208-e0e47529a953"
-version = "0.7.10"
+version = "0.7.11"
 
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SortingNetworks = "07e3d4f1-5dc2-5a3a-9c19-e1965b76eff9"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [weakdeps]
@@ -24,6 +26,7 @@ OffsetArrays = "1"
 Plots = "1"
 Random = "1"
 RecipesBase = "1"
+SortingNetworks = "0.2"
 Statistics = "1"
 Test = "1"
 Unitful = "1"
diff --git a/src/interval.jl b/src/interval.jl
index 101702b..4886ba4 100644
--- a/src/interval.jl
+++ b/src/interval.jl
@@ -161,11 +161,18 @@ show(io::IO, I::OpenInterval) = print(io, leftendpoint(I), " .. ", rightendpoint
 show(io::IO, I::Interval{:open,:closed}) = print(io, leftendpoint(I), " .. ", rightendpoint(I), " (open-closed)")
 show(io::IO, I::Interval{:closed,:open}) = print(io, leftendpoint(I), " .. ", rightendpoint(I), " (closed-open)")
 
+leftendpointtype(::TypedEndpointsInterval{L,R}) where {L,R} = L
+rightendpointtype(::TypedEndpointsInterval{L,R}) where {L,R} = R
+
 # The following are not typestable for mixed endpoint types
 _left_intersect_type(::Type{Val{:open}}, ::Type{Val{L2}}, a1, a2) where L2 = a1 < a2 ? (a2,L2) : (a1,:open)
 _left_intersect_type(::Type{Val{:closed}}, ::Type{Val{L2}}, a1, a2) where L2 = a1 ≤ a2 ? (a2,L2) : (a1,:closed)
 _right_intersect_type(::Type{Val{:open}}, ::Type{Val{R2}}, b1, b2) where R2 = b1 > b2 ? (b2,R2) : (b1,:open)
 _right_intersect_type(::Type{Val{:closed}}, ::Type{Val{R2}}, b1, b2) where R2 = b1 ≥ b2 ? (b2,R2) : (b1,:closed)
+_left_union_type(::Type{Val{:open}}, ::Type{Val{L2}}, a1, a2) where L2 = a1 < a2 ? (a1,:open) : (a2,L2)
+_left_union_type(::Type{Val{:closed}}, ::Type{Val{L2}}, a1, a2) where L2 = a1 ≤ a2 ? (a1,:closed) : (a2,L2)
+_right_union_type(::Type{Val{:open}}, ::Type{Val{R2}}, b1, b2) where R2 = b1 > b2 ? (b1,:open) : (b2,R2)
+_right_union_type(::Type{Val{:closed}}, ::Type{Val{R2}}, b1, b2) where R2 = b1 ≥ b2 ? (b1,:closed) : (b2,R2)
 
 function intersect(d1::TypedEndpointsInterval{L1,R1}, d2::TypedEndpointsInterval{L2,R2}) where {L1,R1,L2,R2}
     a1, b1 = endpoints(d1); a2, b2 = endpoints(d2)
@@ -181,37 +188,12 @@ end
 
 intersect(d1::AbstractInterval, d2::AbstractInterval) = intersect(Interval(d1), Interval(d2))
 
-union(I::TypedEndpointsInterval...) = unioninterval(intervalunion!([I...]))
-unioninterval(x) = isone(length(x)) ? x[1] : throw(ArgumentError("Cannot construct union of disjoint sets."))
-function intervalunion!(I::Vector{<:TypedEndpointsInterval})
-    sort!(I, by = leftendpoint)
-
-    # filter out empty sets
-    k = 1
-    while k <= length(I)
-        if isempty(I[k])
-            popat!(I,k)
-        else
-            k += 1
-        end
-    end
+include("unionalgorithms.jl")
 
-    # merge intervals
-    k = 2
-    while k <= length(I)
-        x = leftendpoint(I[k])
-        y = rightendpoint(I[k-1])
-        if x > y
-            k += 1
-        elseif x < y || x ∈ I[k] || x ∈ I[k-1]
-            I[k-1] = _union(I[k-1],I[k])
-            popat!(I,k)
-        else
-            k += 1
-        end
-    end
-    I
-end
+union(d::TypedEndpointsInterval) = d # 1 interval
+union(d1::TypedEndpointsInterval, d2::TypedEndpointsInterval) = union2(d1, d2) # 2 intervals
+Base.@nexprs(23,N -> union(I::Vararg{TypedEndpointsInterval,N+2}) = tupleunion(swapsort(I))) # 3 to 25 intervals
+union(I::TypedEndpointsInterval...) = tupleunion(sort(SVector(I); lt = leftof)) # ≥26 intervals
 
 # these assume overlap
 function _union(A::TypedEndpointsInterval{L,R}, B::TypedEndpointsInterval{L,R}) where {L,R}
diff --git a/src/unionalgorithms.jl b/src/unionalgorithms.jl
new file mode 100644
index 0000000..b9b389d
--- /dev/null
+++ b/src/unionalgorithms.jl
@@ -0,0 +1,70 @@
+using SortingNetworks # not a permanent solution as the package is not maintained. However it's really fast
+import StaticArrays: SVector
+
+"""
+    leftof(I1::TypedEndpointsInterval, I2::TypedEndpointsInterval)
+
+Returns if `I1` has a part to the left of `I2`.
+"""
+function leftof(I1::TypedEndpointsInterval{L1,R1}, I2::TypedEndpointsInterval{L2,R2}) where {L1,R1,L2,R2}
+    if leftendpoint(I1) < leftendpoint(I2)
+        true
+    elseif leftendpoint(I1) > leftendpoint(I2)
+        false
+    elseif L1 == :closed && L2 == :open
+        true
+    else
+        false
+    end
+end
+
+"""
+    canunion(d1, d2)
+
+Returns if `d1 ∪ d2` is a single interval.
+"""
+@inline canunion(d1, d2) = any(∈(d1), endpoints(d2)) || any(∈(d2), endpoints(d1))
+
+function tupleunion(I::Tuple{Vararg{TypedEndpointsInterval,N}}) where N
+    T = promote_type(map(eltype, I))
+    next = iterate(I)
+    while !isnothing(next)
+        (item, state) = next
+        # find the first non-empty interval
+        if isempty(item)
+            next = iterate(I, state)
+            continue
+        end
+        L = leftendpointtype(item)
+        R = rightendpointtype(item)
+        l = leftendpoint(item)
+        r = rightendpoint(item)
+        next = iterate(I, state)
+        while !isnothing(next)
+            (item, state) = next
+            if isempty(item)
+            elseif leftendpoint(item) > r
+                throw(ArgumentError("IntervalSets doesn't support union of disjoint intervals, while the interval $r..$(leftendpoint(item)) is not covered. Try using DomainSets.UnionDomain for disjoint intervals or ∪(a,b,c...) if the intervals are not sorted."))
+            elseif R==:open && leftendpoint(item)==r && leftendpointtype(item)==:open
+                throw(ArgumentError("IntervalSets doesn't support union of disjoint intervals, while the point $r is not covered. Try using DomainSets.UnionDomain for disjoint intervals or ∪(a,b,c...) if the intervals are not sorted."))
+            else
+                (r,R) = _right_union_type(Val{R}, Val{rightendpointtype(item)}, r, rightendpoint(item))
+            end
+            next = iterate(I, state)
+        end
+        return Interval{L,R,T}(l,r)
+    end
+    return one(T)..zero(T) # can't find the first non-empty interval. return an empty interval.
+end
+# workaround since SortingNetworks.jl doesn't support `lt` keyword argument
+Base.minmax(I1::TypedEndpointsInterval, I2::TypedEndpointsInterval) = leftof(I1,I2) ? (I1,I2) : (I2,I1)
+
+# good old union
+function union2(d1::TypedEndpointsInterval{L1,R1,T1}, d2::TypedEndpointsInterval{L2,R2,T2}) where {L1,R1,T1,L2,R2,T2}
+    T = promote_type(T1,T2)
+    isempty(d1) && return Interval{L2,R2,T}(d2)
+    isempty(d2) && return Interval{L1,R1,T}(d1)
+    any(∈(d1), endpoints(d2)) || any(∈(d2), endpoints(d1)) ||
+        throw(ArgumentError("Cannot construct union of disjoint sets."))
+    _union(d1, d2)
+end

From 638b8e820eb2659c9569484ee18dc9464028974e Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 19 Feb 2024 19:47:57 +0000
Subject: [PATCH 04/17] Update Project.toml

---
 Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Project.toml b/Project.toml
index 882f04a..b45a58f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,6 +27,7 @@ Plots = "1"
 Random = "1"
 RecipesBase = "1"
 SortingNetworks = "0.2"
+StaticArrays = "1"
 Statistics = "1"
 Test = "1"
 Unitful = "1"

From 66d885cd725d7d60ca8491d9920dd679e33393c8 Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 19 Feb 2024 20:43:13 +0000
Subject: [PATCH 05/17] add tests

---
 src/interval.jl        |  4 ++--
 src/unionalgorithms.jl | 12 ++++++------
 test/setoperations.jl  |  6 ++++++
 3 files changed, 14 insertions(+), 8 deletions(-)

diff --git a/src/interval.jl b/src/interval.jl
index 4886ba4..d282fb4 100644
--- a/src/interval.jl
+++ b/src/interval.jl
@@ -192,8 +192,8 @@ include("unionalgorithms.jl")
 
 union(d::TypedEndpointsInterval) = d # 1 interval
 union(d1::TypedEndpointsInterval, d2::TypedEndpointsInterval) = union2(d1, d2) # 2 intervals
-Base.@nexprs(23,N -> union(I::Vararg{TypedEndpointsInterval,N+2}) = tupleunion(swapsort(I))) # 3 to 25 intervals
-union(I::TypedEndpointsInterval...) = tupleunion(sort(SVector(I); lt = leftof)) # ≥26 intervals
+Base.@nexprs(23,N -> union(I::Vararg{TypedEndpointsInterval,N+2}) = iterunion(swapsort(I))) # 3 to 25 intervals
+union(I::TypedEndpointsInterval...) = iterunion(sort!(SVector(I); lt = leftof)) # ≥26 intervals
 
 # these assume overlap
 function _union(A::TypedEndpointsInterval{L,R}, B::TypedEndpointsInterval{L,R}) where {L,R}
diff --git a/src/unionalgorithms.jl b/src/unionalgorithms.jl
index b9b389d..74078b5 100644
--- a/src/unionalgorithms.jl
+++ b/src/unionalgorithms.jl
@@ -25,21 +25,21 @@ Returns if `d1 ∪ d2` is a single interval.
 """
 @inline canunion(d1, d2) = any(∈(d1), endpoints(d2)) || any(∈(d2), endpoints(d1))
 
-function tupleunion(I::Tuple{Vararg{TypedEndpointsInterval,N}}) where N
-    T = promote_type(map(eltype, I))
-    next = iterate(I)
+function iterunion(iter)
+    T = promote_type(map(eltype, iter)...)
+    next = iterate(iter)
     while !isnothing(next)
         (item, state) = next
         # find the first non-empty interval
         if isempty(item)
-            next = iterate(I, state)
+            next = iterate(iter, state)
             continue
         end
         L = leftendpointtype(item)
         R = rightendpointtype(item)
         l = leftendpoint(item)
         r = rightendpoint(item)
-        next = iterate(I, state)
+        next = iterate(iter, state)
         while !isnothing(next)
             (item, state) = next
             if isempty(item)
@@ -50,7 +50,7 @@ function tupleunion(I::Tuple{Vararg{TypedEndpointsInterval,N}}) where N
             else
                 (r,R) = _right_union_type(Val{R}, Val{rightendpointtype(item)}, r, rightendpoint(item))
             end
-            next = iterate(I, state)
+            next = iterate(iter, state)
         end
         return Interval{L,R,T}(l,r)
     end
diff --git a/test/setoperations.jl b/test/setoperations.jl
index 551300c..29e7347 100644
--- a/test/setoperations.jl
+++ b/test/setoperations.jl
@@ -230,6 +230,12 @@
     # - different interval types
     @test (1..2) ∩ OpenInterval(0.5, 1.5) ≡ Interval{:closed, :open}(1, 1.5)
     @test (1..2) ∪ OpenInterval(0.5, 1.5) ≡ Interval{:open, :closed}(0.5, 2)
+
+    # union of multiple intervals
+    intervals = [i1, i2, i3, i4, i5, i_empty]
+    for _ in 1:10
+        @test ∪(shuffle!(intervals)...) == 0..3
+    end
 end
 
 @testset "in" begin

From 22803b98a818245421b8af58103bceed108eaba9 Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 19 Feb 2024 20:59:58 +0000
Subject: [PATCH 06/17] add tests

---
 test/setoperations.jl | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/test/setoperations.jl b/test/setoperations.jl
index 29e7347..38c955f 100644
--- a/test/setoperations.jl
+++ b/test/setoperations.jl
@@ -236,6 +236,14 @@
     for _ in 1:10
         @test ∪(shuffle!(intervals)...) == 0..3
     end
+    intervals = [i1, i2, i4, i5, i_empty]
+    for _ in 1:10
+        @test_throws ArgumentError ∪(shuffle!(intervals)...)
+    end
+    intervals = [OpenInterval(n-0.5,n+0.5) for n in 1:50]
+    for _ in 1:10
+        @test ∪(shuffle!(intervals)...) == 1..0
+    end
 end
 
 @testset "in" begin

From e53a3d9e0f0bd733c8aeb4909f7abf72f8a0dd94 Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 19 Feb 2024 22:37:35 +0000
Subject: [PATCH 07/17] Update setoperations.jl

---
 test/setoperations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/setoperations.jl b/test/setoperations.jl
index 38c955f..a914e55 100644
--- a/test/setoperations.jl
+++ b/test/setoperations.jl
@@ -238,7 +238,7 @@
     end
     intervals = [i1, i2, i4, i5, i_empty]
     for _ in 1:10
-        @test_throws ArgumentError ∪(shuffle!(intervals)...)
+        @test_throws ArgumentError union(shuffle!(intervals)...)
     end
     intervals = [OpenInterval(n-0.5,n+0.5) for n in 1:50]
     for _ in 1:10

From 4aaf4d093beb69b7831542cfb01cbf21c3de156a Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 19 Feb 2024 23:08:49 +0000
Subject: [PATCH 08/17] Update interval.jl

---
 src/interval.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/interval.jl b/src/interval.jl
index d282fb4..aae1841 100644
--- a/src/interval.jl
+++ b/src/interval.jl
@@ -193,7 +193,7 @@ include("unionalgorithms.jl")
 union(d::TypedEndpointsInterval) = d # 1 interval
 union(d1::TypedEndpointsInterval, d2::TypedEndpointsInterval) = union2(d1, d2) # 2 intervals
 Base.@nexprs(23,N -> union(I::Vararg{TypedEndpointsInterval,N+2}) = iterunion(swapsort(I))) # 3 to 25 intervals
-union(I::TypedEndpointsInterval...) = iterunion(sort!(SVector(I); lt = leftof)) # ≥26 intervals
+union(I::TypedEndpointsInterval...) = iterunion(sort(SVector(I); lt = leftof)) # ≥26 intervals
 
 # these assume overlap
 function _union(A::TypedEndpointsInterval{L,R}, B::TypedEndpointsInterval{L,R}) where {L,R}

From c26b46672a58294bc07514667fbb092eb3b70e3a Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Mon, 19 Feb 2024 23:15:31 +0000
Subject: [PATCH 09/17] Update setoperations.jl

---
 test/setoperations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/setoperations.jl b/test/setoperations.jl
index a914e55..ff228a3 100644
--- a/test/setoperations.jl
+++ b/test/setoperations.jl
@@ -240,7 +240,7 @@
     for _ in 1:10
         @test_throws ArgumentError union(shuffle!(intervals)...)
     end
-    intervals = [OpenInterval(n-0.5,n+0.5) for n in 1:50]
+    intervals = [OpenInterval(n+0.5,n-0.5) for n in 1:50]
     for _ in 1:10
         @test ∪(shuffle!(intervals)...) == 1..0
     end

From 2d30255bbb03de5e6bd87074ed3dfefdeec589db Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Tue, 20 Feb 2024 17:10:08 +0000
Subject: [PATCH 10/17] change to TupleTools and add tests

---
 Project.toml           | 10 +++++---
 src/interval.jl        |  4 ++--
 src/unionalgorithms.jl | 11 ++++-----
 test/runtests.jl       |  2 ++
 test/setoperations.jl  | 53 ++++++++++++++++++++++++++++++++++++++----
 5 files changed, 64 insertions(+), 16 deletions(-)

diff --git a/Project.toml b/Project.toml
index b45a58f..7292896 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,9 +5,9 @@ version = "0.7.11"
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
-SortingNetworks = "07e3d4f1-5dc2-5a3a-9c19-e1965b76eff9"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 
 [weakdeps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -22,27 +22,31 @@ IntervalSetsStatisticsExt = "Statistics"
 [compat]
 Aqua = "0.8"
 Dates = "1"
+DomainSets = "0.7"
 OffsetArrays = "1"
 Plots = "1"
 Random = "1"
 RecipesBase = "1"
-SortingNetworks = "0.2"
+Scanf = "0.5"
 StaticArrays = "1"
 Statistics = "1"
 Test = "1"
+TupleTools = "1"
 Unitful = "1"
 julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
+DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+Scanf = "6ef1bc8b-493b-44e1-8d40-549aa65c4b41"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["Aqua", "Dates", "Test", "Plots", "Random", "RecipesBase", "OffsetArrays", "Statistics", "Unitful"]
+test = ["Aqua", "Dates", "DomainSets", "Test", "Plots", "Random", "RecipesBase", "OffsetArrays", "Scanf", "Statistics", "Unitful"]
diff --git a/src/interval.jl b/src/interval.jl
index aae1841..6d9adef 100644
--- a/src/interval.jl
+++ b/src/interval.jl
@@ -192,8 +192,8 @@ include("unionalgorithms.jl")
 
 union(d::TypedEndpointsInterval) = d # 1 interval
 union(d1::TypedEndpointsInterval, d2::TypedEndpointsInterval) = union2(d1, d2) # 2 intervals
-Base.@nexprs(23,N -> union(I::Vararg{TypedEndpointsInterval,N+2}) = iterunion(swapsort(I))) # 3 to 25 intervals
-union(I::TypedEndpointsInterval...) = iterunion(sort(SVector(I); lt = leftof)) # ≥26 intervals
+Base.@nexprs(4,N -> union(I::Vararg{TypedEndpointsInterval,N+2}) = iterunion(TupleTools.sort(I; lt = leftof))) # 3 to 6 intervals
+union(I::TypedEndpointsInterval...) = iterunion(sort(SVector(I); lt = leftof).data) # ≥7 intervals
 
 # these assume overlap
 function _union(A::TypedEndpointsInterval{L,R}, B::TypedEndpointsInterval{L,R}) where {L,R}
diff --git a/src/unionalgorithms.jl b/src/unionalgorithms.jl
index 74078b5..f5d0b23 100644
--- a/src/unionalgorithms.jl
+++ b/src/unionalgorithms.jl
@@ -1,4 +1,4 @@
-using SortingNetworks # not a permanent solution as the package is not maintained. However it's really fast
+import TupleTools
 import StaticArrays: SVector
 
 """
@@ -44,7 +44,7 @@ function iterunion(iter)
             (item, state) = next
             if isempty(item)
             elseif leftendpoint(item) > r
-                throw(ArgumentError("IntervalSets doesn't support union of disjoint intervals, while the interval $r..$(leftendpoint(item)) is not covered. Try using DomainSets.UnionDomain for disjoint intervals or ∪(a,b,c...) if the intervals are not sorted."))
+                throw(ArgumentError("IntervalSets doesn't support union of disjoint intervals, while the interval $r..$(leftendpoint(item)) (open) is not covered. Try using DomainSets.UnionDomain for disjoint intervals or ∪(a,b,c...) if the intervals are not sorted."))
             elseif R==:open && leftendpoint(item)==r && leftendpointtype(item)==:open
                 throw(ArgumentError("IntervalSets doesn't support union of disjoint intervals, while the point $r is not covered. Try using DomainSets.UnionDomain for disjoint intervals or ∪(a,b,c...) if the intervals are not sorted."))
             else
@@ -56,15 +56,12 @@ function iterunion(iter)
     end
     return one(T)..zero(T) # can't find the first non-empty interval. return an empty interval.
 end
-# workaround since SortingNetworks.jl doesn't support `lt` keyword argument
-Base.minmax(I1::TypedEndpointsInterval, I2::TypedEndpointsInterval) = leftof(I1,I2) ? (I1,I2) : (I2,I1)
 
 # good old union
 function union2(d1::TypedEndpointsInterval{L1,R1,T1}, d2::TypedEndpointsInterval{L2,R2,T2}) where {L1,R1,T1,L2,R2,T2}
     T = promote_type(T1,T2)
     isempty(d1) && return Interval{L2,R2,T}(d2)
     isempty(d2) && return Interval{L1,R1,T}(d1)
-    any(∈(d1), endpoints(d2)) || any(∈(d2), endpoints(d1)) ||
-        throw(ArgumentError("Cannot construct union of disjoint sets."))
-    _union(d1, d2)
+    canunion(d1, d2) && return _union(d1, d2)
+    throw(ArgumentError("Cannot construct union of disjoint sets."))
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index d9dc0ea..e620fe0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,6 +6,8 @@ import Statistics: mean
 using Random
 using Unitful
 using Plots
+using Scanf
+using DomainSets
 
 import IntervalSets: Domain, endpoints, closedendpoints, TypedEndpointsInterval
 
diff --git a/test/setoperations.jl b/test/setoperations.jl
index ff228a3..3bcdd83 100644
--- a/test/setoperations.jl
+++ b/test/setoperations.jl
@@ -231,7 +231,6 @@
     @test (1..2) ∩ OpenInterval(0.5, 1.5) ≡ Interval{:closed, :open}(1, 1.5)
     @test (1..2) ∪ OpenInterval(0.5, 1.5) ≡ Interval{:open, :closed}(0.5, 2)
 
-    # union of multiple intervals
     intervals = [i1, i2, i3, i4, i5, i_empty]
     for _ in 1:10
         @test ∪(shuffle!(intervals)...) == 0..3
@@ -240,9 +239,55 @@
     for _ in 1:10
         @test_throws ArgumentError union(shuffle!(intervals)...)
     end
-    intervals = [OpenInterval(n+0.5,n-0.5) for n in 1:50]
-    for _ in 1:10
-        @test ∪(shuffle!(intervals)...) == 1..0
+end
+
+randinterval(s) = Interval{rand([:closed,:open]),rand([:closed,:open])}(rand(s), rand(s))
+function test_multipleunion(intervals)
+    if all(isempty, intervals)
+        @test ∪(intervals...) == 1 .. 0
+    else
+        u = nothing
+        try
+            u = ∪(intervals...)
+        catch e
+            @test e isa ArgumentError
+            s = e.msg
+            ind = findfirst("while the ", s)[end] + 1
+            if startswith(s[ind:end], "interval")
+                u = OpenInterval(@scanf(s[ind+9:end], "%d..%d", Int, Int)[2:3]...)
+                @test all(v -> isempty(u ∩ v), intervals)
+            elseif startswith(s[ind:end], "point")
+                x = @scanf(s[ind+6:end], "%d", Int)[2]
+                @test all(v -> x ∉ v, intervals)
+            else
+                error("have you touched the error message?")
+            end
+            return
+        end
+        @test all(v -> v ⊆ u, intervals)
+        # due to an issue from DomainSets.jl, setdiff is currently unreliable. See github.com/JuliaApproximation/DomainSets.jl/issues/151
+        # as a result, the correctness of interval union is not thoroughly tested.
+        #= v = u
+        for i in intervals
+            try
+                v = setdiff(v, i)
+            catch
+                println("setdiff($v, $i) was not successful. The union was $u and the components are $intervals.")
+                @test false
+            end
+        end
+        @test isempty(v) =#
+    end
+end
+
+@testset "general union" begin
+    for _ in 1:500
+        intervals = [randinterval(0:10) for _ in 1:5]
+        test_multipleunion(intervals)
+    end
+    for _ in 1:50
+        intervals = [randinterval(0:100) for _ in 1:100]
+        test_multipleunion(intervals)
     end
 end
 

From 383dd15fde61e5918c292111a4ca41639e929386 Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Thu, 22 Feb 2024 09:55:20 +0000
Subject: [PATCH 11/17] explain `canunion` better

---
 src/unionalgorithms.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/unionalgorithms.jl b/src/unionalgorithms.jl
index f5d0b23..5063781 100644
--- a/src/unionalgorithms.jl
+++ b/src/unionalgorithms.jl
@@ -21,7 +21,7 @@ end
 """
     canunion(d1, d2)
 
-Returns if `d1 ∪ d2` is a single interval.
+Returns if `d1 ∪ d2` is a single interval. `d1` and `d2` have to be non-empty. Note that `canunion` is not always `!isdisjoint`. For example, ``[0,1)`` and ``[1,2]`` are disjoint, but the union of them is a single interval.
 """
 @inline canunion(d1, d2) = any(∈(d1), endpoints(d2)) || any(∈(d2), endpoints(d1))
 

From f33a621978d61be503af2a616df37a909936624d Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Thu, 22 Feb 2024 09:56:06 +0000
Subject: [PATCH 12/17] move to StaticArraysCore

---
 Project.toml           | 4 ++--
 src/unionalgorithms.jl | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/Project.toml b/Project.toml
index 7292896..46b7902 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,7 +5,7 @@ version = "0.7.11"
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 
@@ -28,7 +28,7 @@ Plots = "1"
 Random = "1"
 RecipesBase = "1"
 Scanf = "0.5"
-StaticArrays = "1"
+StaticArraysCore = "1"
 Statistics = "1"
 Test = "1"
 TupleTools = "1"
diff --git a/src/unionalgorithms.jl b/src/unionalgorithms.jl
index 5063781..e458ac8 100644
--- a/src/unionalgorithms.jl
+++ b/src/unionalgorithms.jl
@@ -1,5 +1,5 @@
 import TupleTools
-import StaticArrays: SVector
+import StaticArraysCore: SVector
 
 """
     leftof(I1::TypedEndpointsInterval, I2::TypedEndpointsInterval)

From 7e5668b0f86da1a8f08117eb8dbff5404481e184 Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Thu, 22 Feb 2024 09:59:24 +0000
Subject: [PATCH 13/17] Update unionalgorithms.jl

---
 src/unionalgorithms.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/unionalgorithms.jl b/src/unionalgorithms.jl
index e458ac8..93b3b67 100644
--- a/src/unionalgorithms.jl
+++ b/src/unionalgorithms.jl
@@ -54,7 +54,7 @@ function iterunion(iter)
         end
         return Interval{L,R,T}(l,r)
     end
-    return one(T)..zero(T) # can't find the first non-empty interval. return an empty interval.
+    return first(iter) # can't find the first non-empty interval. return the first one.
 end
 
 # good old union

From 060e1fe911f4bc290f9c6998aa514b4fdef0a0cc Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Thu, 22 Feb 2024 10:00:02 +0000
Subject: [PATCH 14/17] Update setoperations.jl

---
 test/setoperations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/setoperations.jl b/test/setoperations.jl
index 3bcdd83..4b99582 100644
--- a/test/setoperations.jl
+++ b/test/setoperations.jl
@@ -244,7 +244,7 @@ end
 randinterval(s) = Interval{rand([:closed,:open]),rand([:closed,:open])}(rand(s), rand(s))
 function test_multipleunion(intervals)
     if all(isempty, intervals)
-        @test ∪(intervals...) == 1 .. 0
+        @test isempty(∪(intervals...))
     else
         u = nothing
         try

From 332bbf244371bfd4448c8a95062251422ec6429f Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Thu, 22 Feb 2024 10:07:29 +0000
Subject: [PATCH 15/17] explain why the tests are removed

---
 test/setoperations.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/test/setoperations.jl b/test/setoperations.jl
index 4b99582..226c565 100644
--- a/test/setoperations.jl
+++ b/test/setoperations.jl
@@ -265,7 +265,9 @@ function test_multipleunion(intervals)
             return
         end
         @test all(v -> v ⊆ u, intervals)
-        # due to an issue from DomainSets.jl, setdiff is currently unreliable. See github.com/JuliaApproximation/DomainSets.jl/issues/151
+        # The following codes rigorously tests the correctness of u. 
+        # However, the `setdiff` is only implemented in `DomainSets.jl`
+        # which introduces piracies. See https://github.com/JuliaMath/IntervalSets.jl/pull/156#discussion_r1497829695
         # as a result, the correctness of interval union is not thoroughly tested.
         #= v = u
         for i in intervals

From 739ee4464208819c128edf10018dbb197158f3ef Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Thu, 22 Feb 2024 10:54:30 +0000
Subject: [PATCH 16/17] remove DomainSets

---
 Project.toml     | 4 +---
 test/runtests.jl | 1 -
 2 files changed, 1 insertion(+), 4 deletions(-)

diff --git a/Project.toml b/Project.toml
index 46b7902..0a1a67c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -22,7 +22,6 @@ IntervalSetsStatisticsExt = "Statistics"
 [compat]
 Aqua = "0.8"
 Dates = "1"
-DomainSets = "0.7"
 OffsetArrays = "1"
 Plots = "1"
 Random = "1"
@@ -38,7 +37,6 @@ julia = "1.6"
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
-DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -49,4 +47,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["Aqua", "Dates", "DomainSets", "Test", "Plots", "Random", "RecipesBase", "OffsetArrays", "Scanf", "Statistics", "Unitful"]
+test = ["Aqua", "Dates", "Test", "Plots", "Random", "RecipesBase", "OffsetArrays", "Scanf", "Statistics", "Unitful"]
diff --git a/test/runtests.jl b/test/runtests.jl
index e620fe0..189a88b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,7 +7,6 @@ using Random
 using Unitful
 using Plots
 using Scanf
-using DomainSets
 
 import IntervalSets: Domain, endpoints, closedendpoints, TypedEndpointsInterval
 

From a099521bf99e7521226157692155bcec0169c807 Mon Sep 17 00:00:00 2001
From: Tianyi Pu <912396513@qq.com>
Date: Thu, 22 Feb 2024 13:03:28 +0000
Subject: [PATCH 17/17] remove StaticArraysCore

---
 Project.toml           |  2 --
 src/interval.jl        | 33 ++-------------------------------
 src/unionalgorithms.jl |  8 +++++++-
 3 files changed, 9 insertions(+), 34 deletions(-)

diff --git a/Project.toml b/Project.toml
index 0a1a67c..d04f82c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,7 +5,6 @@ version = "0.7.11"
 [deps]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
-StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6"
 
@@ -27,7 +26,6 @@ Plots = "1"
 Random = "1"
 RecipesBase = "1"
 Scanf = "0.5"
-StaticArraysCore = "1"
 Statistics = "1"
 Test = "1"
 TupleTools = "1"
diff --git a/src/interval.jl b/src/interval.jl
index 6d9adef..2cadebf 100644
--- a/src/interval.jl
+++ b/src/interval.jl
@@ -192,37 +192,8 @@ include("unionalgorithms.jl")
 
 union(d::TypedEndpointsInterval) = d # 1 interval
 union(d1::TypedEndpointsInterval, d2::TypedEndpointsInterval) = union2(d1, d2) # 2 intervals
-Base.@nexprs(4,N -> union(I::Vararg{TypedEndpointsInterval,N+2}) = iterunion(TupleTools.sort(I; lt = leftof))) # 3 to 6 intervals
-union(I::TypedEndpointsInterval...) = iterunion(sort(SVector(I); lt = leftof).data) # ≥7 intervals
-
-# these assume overlap
-function _union(A::TypedEndpointsInterval{L,R}, B::TypedEndpointsInterval{L,R}) where {L,R}
-    left = min(leftendpoint(A), leftendpoint(B))
-    right = max(rightendpoint(A), rightendpoint(B))
-    Interval{L,R}(left, right)
-end
-
-# this is not typestable
-function _union(A::TypedEndpointsInterval{L1,R1}, B::TypedEndpointsInterval{L2,R2}) where {L1,R1,L2,R2}
-    if leftendpoint(A) == leftendpoint(B)
-        L = L1 == :closed ? :closed : L2
-    elseif leftendpoint(A) < leftendpoint(B)
-        L = L1
-    else
-        L = L2
-    end
-    if rightendpoint(A) == rightendpoint(B)
-        R = R1 == :closed ? :closed : R2
-    elseif rightendpoint(A) > rightendpoint(B)
-        R = R1
-    else
-        R = R2
-    end
-    left = min(leftendpoint(A), leftendpoint(B))
-    right = max(rightendpoint(A), rightendpoint(B))
-
-    Interval{L,R}(left, right)
-end
+Base.@nexprs(18,N -> union(I::Vararg{TypedEndpointsInterval,N+2}) = iterunion(TupleTools.sort(I; lt = leftof))) # 3 to 20 intervals
+union(I::TypedEndpointsInterval...) = iterunion(sort!(collect(I); lt = leftof)) # ≥21 intervals
 
 ClosedInterval{T}(i::AbstractUnitRange{I}) where {T,I<:Integer} = ClosedInterval{T}(minimum(i), maximum(i))
 ClosedInterval(i::AbstractUnitRange{I}) where {I<:Integer} = ClosedInterval{I}(minimum(i), maximum(i))
diff --git a/src/unionalgorithms.jl b/src/unionalgorithms.jl
index 93b3b67..8994653 100644
--- a/src/unionalgorithms.jl
+++ b/src/unionalgorithms.jl
@@ -1,5 +1,4 @@
 import TupleTools
-import StaticArraysCore: SVector
 
 """
     leftof(I1::TypedEndpointsInterval, I2::TypedEndpointsInterval)
@@ -65,3 +64,10 @@ function union2(d1::TypedEndpointsInterval{L1,R1,T1}, d2::TypedEndpointsInterval
     canunion(d1, d2) && return _union(d1, d2)
     throw(ArgumentError("Cannot construct union of disjoint sets."))
 end
+
+# this is not typestable
+function _union(A::TypedEndpointsInterval{L1,R1}, B::TypedEndpointsInterval{L2,R2}) where {L1,R1,L2,R2}
+    l, L = _left_union_type(Val{L1}, Val{L2}, leftendpoint(A), leftendpoint(B))
+    r, R = _right_union_type(Val{R1}, Val{R2}, rightendpoint(A), rightendpoint(B))
+    Interval{L,R}(l, r)
+end
\ No newline at end of file