add intersect for CartesianIndices (#36643)

NEWS.md

@@ -67,6 +67,7 @@ Standard library changes
 * `Libdl` has been moved to `Base.Libc.Libdl`, however it is still accessible as an stdlib ([#35628]).
 * `first` and `last` functions now accept an integer as second argument to get that many
   leading or trailing elements of any iterable ([#34868]).
+* `intersect` on `CartesianIndices` now returns `CartesianIndices` instead of `Vector{<:CartesianIndex}` ([#36643]).
 
 #### LinearAlgebra
 * New method `LinearAlgebra.issuccess(::CholeskyPivoted)` for checking whether pivoted Cholesky factorization was successful ([#36002]).

base/multidimensional.jl

@@ -479,6 +479,10 @@ module IteratorsMD
 
     Base.LinearIndices(inds::CartesianIndices{N,R}) where {N,R} = LinearIndices{N,R}(inds.indices)
 
+    # array operations
+    Base.intersect(a::CartesianIndices{N}, b::CartesianIndices{N}) where N =
+        CartesianIndices(intersect.(a.indices, b.indices))
+
     # Views of reshaped CartesianIndices are used for partitions — ensure these are fast
     const CartesianPartition{T<:CartesianIndex, P<:CartesianIndices, R<:ReshapedArray{T,1,P}} = SubArray{T,1,R,Tuple{UnitRange{Int}},false}
     eltype(::Type{PartitionIterator{T}}) where {T<:ReshapedArrayLF} = SubArray{eltype(T), 1, T, Tuple{UnitRange{Int}}, true}

test/cartesian.jl

@@ -67,4 +67,9 @@ end
     @test iterate(I, CartesianIndex(4, typemax(Int)))    === nothing
 end
 
+@testset "CartesianIndices operations" begin
+    I = CartesianIndices((1:3, 4:6))
+    J = CartesianIndices((2:4, 3:5))
 
+    @test @inferred(intersect(I, J)) == CartesianIndices((2:3, 4:5))
+end