diff --git a/test/PolyhedralGeometry/polyhedron.jl b/test/PolyhedralGeometry/polyhedron.jl
index 8b6d15c0d1da..1445276d7405 100644
--- a/test/PolyhedralGeometry/polyhedron.jl
+++ b/test/PolyhedralGeometry/polyhedron.jl
@@ -37,9 +37,17 @@
   @test nrows(Oscar.pm_object(pf3).INEQUALITIES) == 4
   @test n_vertices(pf3) == 1
 
+  function _check_im_perm_rows(inc::IncidenceMatrix, o)
+    oinc = IncidenceMatrix(o)
+    nr = nrows(inc)
+    nr == nrows(oinc) &&
+      issetequal(Polymake.row.(Ref(inc), 1:nr),
+                 Polymake.row.(Ref(oinc), 1:nr))
+  end
+
   @testset "core functionality" begin
     @test matrix(f, rays(Q1)) * v == T[f(2)]
-    @test matrix(f, vertices(Q1)) * v == T[f(1), f(0), f(1)]
+    @test issetequal(matrix(f, vertices(Q1)) * v, T[f(1), f(0), f(1)])
     @test issubset(Q0, Q1)
     @test !issubset(Q1, Q0)
     @test [1, 0] in Q0
@@ -47,21 +55,20 @@
     @test n_vertices(Q0) == 3
     @test n_vertices.(faces(Q0, 1)) == [2, 2, 2]
       @test lattice_points(Q0) isa SubObjectIterator{PointVector{ZZRingElem}}
-      @test point_matrix(lattice_points(Q0)) == matrix(ZZ, [0 0; 0 1; 1 0])
-      @test matrix(ZZ, lattice_points(Q0)) == matrix(ZZ, [0 0; 0 1; 1 0])
+      @test point_matrix(lattice_points(Q0)) isa MatElem{ZZRingElem}
+      @test matrix(ZZ, lattice_points(Q0)) isa MatElem{ZZRingElem}
       @test length(lattice_points(Q0)) == 3
-      @test lattice_points(Q0) == [[0, 0], [0, 1], [1, 0]]
+      @test issetequal(lattice_points(Q0), point_vector.(Ref(ZZ), [[0, 0], [0, 1], [1, 0]]))
       @test interior_lattice_points(square) isa SubObjectIterator{PointVector{ZZRingElem}}
       @test point_matrix(interior_lattice_points(square)) == matrix(ZZ, [0 0])
       @test matrix(ZZ, interior_lattice_points(square)) == matrix(ZZ, [0 0])
       @test length(interior_lattice_points(square)) == 1
       @test interior_lattice_points(square) == [[0, 0]]
       @test boundary_lattice_points(square) isa SubObjectIterator{PointVector{ZZRingElem}}
-      @test point_matrix(boundary_lattice_points(square)) ==
-        matrix(ZZ, [-1 -1; -1 0; -1 1; 0 -1; 0 1; 1 -1; 1 0; 1 1])
+      @test point_matrix(boundary_lattice_points(square)) isa MatElem{ZZRingElem}
       @test length(boundary_lattice_points(square)) == 8
-      @test boundary_lattice_points(square) ==
-        [[-1, -1], [-1, 0], [-1, 1], [0, -1], [0, 1], [1, -1], [1, 0], [1, 1]]
+      @test issetequal(boundary_lattice_points(square),
+                           point_vector.(Ref(ZZ), [[-1, -1], [-1, 0], [-1, 1], [0, -1], [0, 1], [1, -1], [1, 0], [1, 1]]))
     if T == QQFieldElem
       @test is_smooth(Q0)
       @test is_normal(Q0)
@@ -104,13 +111,8 @@
     @test rays(RayVector, Pos) isa SubObjectIterator{RayVector{T}}
     @test rays(Pos) isa SubObjectIterator{RayVector{T}}
     @test length(rays(Pos)) == 3
-    if T == QQFieldElem
-      @test vector_matrix(rays(Pos)) == matrix(QQ, [1 0 0; 0 1 0; 0 0 1])
-      @test rays(Pos) == [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
-    else
-      @test vector_matrix(rays(Pos)) == matrix(f, [0 1 0; 1 0 0; 0 0 1])
-      @test rays(Pos) == [[0, 1, 0], [1, 0, 0], [0, 0, 1]]
-    end
+    @test vector_matrix(rays(Pos)) isa MatElem{T}
+    @test issetequal(rays(Pos), ray_vector.(Ref(f), [[1, 0, 0], [0, 1, 0], [0, 0, 1]]))
     @test lineality_space(L) isa SubObjectIterator{RayVector{T}}
     @test generator_matrix(lineality_space(L)) == matrix(f, [0 0 1])
     if T == QQFieldElem
@@ -120,49 +122,46 @@
     @test lineality_space(L) == [[0, 0, 1]]
     @test faces(square, 1) isa SubObjectIterator{Polyhedron{T}}
     @test length(faces(square, 1)) == 4
-    @test faces(square, 1) ==
-      convex_hull.([f], [[-1 -1; -1 1], [1 -1; 1 1], [-1 -1; 1 -1], [-1 1; 1 1]])
-    @test vertex_indices(faces(square, 1)) ==
-      IncidenceMatrix([[1, 3], [2, 4], [1, 2], [3, 4]])
+    @test issetequal(faces(square, 1),
+                     convex_hull.(Ref(f), [[-1 -1; -1 1], [1 -1; 1 1], [-1 -1; 1 -1], [-1 1; 1 1]]))
+    @test _check_im_perm_rows(vertex_indices(faces(square, 1)),
+                              [[1, 3], [2, 4], [1, 2], [3, 4]])
     @test ray_indices(faces(square, 1)) == IncidenceMatrix(4, 0)
-    @test vertex_and_ray_indices(faces(square, 1)) ==
-      IncidenceMatrix([[1, 3], [2, 4], [1, 2], [3, 4]])
-    @test IncidenceMatrix(faces(square, 1)) ==
-      IncidenceMatrix([[1, 3], [2, 4], [1, 2], [3, 4]])
-    @test faces(IncidenceMatrix, square, 1) ==
-      IncidenceMatrix([[1, 3], [2, 4], [1, 2], [3, 4]])
-    @test facet_indices(vertices(square)) ==
-      IncidenceMatrix([[1, 3], [2, 3], [1, 4], [2, 4]])
-    @test IncidenceMatrix(vertices(square)) ==
-      IncidenceMatrix([[1, 3], [2, 3], [1, 4], [2, 4]])
-    @test vertices(IncidenceMatrix, square) ==
-      IncidenceMatrix([[1, 3], [2, 3], [1, 4], [2, 4]])
+    @test _check_im_perm_rows(vertex_and_ray_indices(faces(square, 1)),
+                              [[2, 4], [1, 3], [1, 2], [3, 4]])
+    @test _check_im_perm_rows(IncidenceMatrix(faces(square, 1)),
+                              [[1, 3], [2, 4], [1, 2], [3, 4]])
+    @test _check_im_perm_rows(faces(IncidenceMatrix, square, 1),
+                              [[1, 3], [2, 4], [1, 2], [3, 4]])
+    @test _check_im_perm_rows(facet_indices(vertices(square)),
+                              [[1, 3], [2, 3], [1, 4], [2, 4]])
+    @test _check_im_perm_rows(IncidenceMatrix(vertices(square)),
+                              [[1, 3], [2, 3], [1, 4], [2, 4]])
+    @test _check_im_perm_rows(vertices(IncidenceMatrix, square),
+                              [[1, 3], [2, 3], [1, 4], [2, 4]])
     @test faces(Pos, 1) isa SubObjectIterator{Polyhedron{T}}
     @test length(faces(Pos, 1)) == 3
-    if T == QQFieldElem
-      @test faces(Pos, 1) == convex_hull.(T, [[0 0 0]], [[1 0 0], [0 1 0], [0 0 1]])
-    else
-      @test faces(Pos, 1) == convex_hull.([f], [[0 0 0]], [[0 1 0], [1 0 0], [0 0 1]])
-    end
-    @test vertex_indices(faces(Pos, 1)) == IncidenceMatrix([[1], [1], [1]])
-    @test ray_indices(faces(Pos, 1)) == IncidenceMatrix([[1], [2], [3]])
-    @test vertex_and_ray_indices(faces(Pos, 1)) == IncidenceMatrix([[1, 4], [2, 4], [3, 4]])
-    @test IncidenceMatrix(faces(Pos, 1)) == IncidenceMatrix([[1, 4], [2, 4], [3, 4]])
-    @test faces(IncidenceMatrix, Pos, 1) == IncidenceMatrix([[1, 4], [2, 4], [3, 4]])
+    @test issetequal(faces(Pos, 1), convex_hull.(Ref(f), [[0 0 0]], [[1 0 0], [0 1 0], [0 0 1]]))
+    @test _check_im_perm_rows(vertex_indices(faces(Pos, 1)), [[1], [1], [1]])
+    @test _check_im_perm_rows(ray_indices(faces(Pos, 1)), [[1], [2], [3]])
+    @test _check_im_perm_rows(vertex_and_ray_indices(faces(Pos, 1)), [[1, 4], [2, 4], [3, 4]])
+    @test _check_im_perm_rows(IncidenceMatrix(faces(Pos, 1)), [[1, 4], [2, 4], [3, 4]])
+    @test _check_im_perm_rows(faces(IncidenceMatrix, Pos, 1), [[1, 4], [2, 4], [3, 4]])
     @test isnothing(faces(Q2, 0))
     v = vertices(minkowski_sum(Q0, square))
     @test length(v) == 5
-    @test v == [f.([2, -1]), f.([2, 1]), f.([-1, -1]), f.([-1, 2]), f.([1, 2])]
-    @test point_matrix(v) == matrix(f, [2 -1; 2 1; -1 -1; -1 2; 1 2])
+    @test issetequal(v, point_vector.(Ref(f), [[2, -1], [2, 1], [-1, -1], [-1, 2], [1, 2]]))
+    @test point_matrix(v) isa MatElem{T}
+    let S = Pair{Matrix{T},T}
+      @test issetequal(facets(S, Pos), S.([f.([-1 0 0]), f.([0 -1 0]), f.([0 0 -1])], [f(0)]))
+    end
+    let S = AffineHalfspace{T}
+      @test issetequal(facets(S, Pos), affine_halfspace.(Ref(f), [[-1 0 0], [0 -1 0], [0 0 -1]], [0]))
+    end
     for S in [AffineHalfspace{T},
-      Pair{Matrix{T},T},
-      Polyhedron{T}]
+              Pair{Matrix{T},T},
+              Polyhedron{T}]
       @test facets(S, Pos) isa SubObjectIterator{S}
-      if S == Pair{Matrix{T},T}
-        @test facets(S, Pos) == S.([f.([-1 0 0]), f.([0 -1 0]), f.([0 0 -1])], [f(0)])
-      elseif S == AffineHalfspace{T}
-        @test facets(S, Pos) == affine_halfspace.([f], [[-1 0 0], [0 -1 0], [0 0 -1]], [0])
-      end
       @test length(facets(S, Pos)) == 3
       @test affine_inequality_matrix(facets(S, Pos)) ==
         matrix(f, [0 -1 0 0; 0 0 -1 0; 0 0 0 -1])
@@ -173,64 +172,31 @@
         @test hmp.A == matrix(f, [-1 0 0; 0 -1 0; 0 0 -1])
         @test hmp.b == f.([0, 0, 0])
       end
-      if T == QQFieldElem
-        @test ray_indices(facets(S, Pos)) == IncidenceMatrix([[2, 3], [1, 3], [1, 2]])
-        @test vertex_and_ray_indices(facets(S, Pos)) ==
-          IncidenceMatrix([[2, 3, 4], [1, 3, 4], [1, 2, 4]])
-        @test IncidenceMatrix(facets(S, Pos)) ==
-          IncidenceMatrix([[2, 3, 4], [1, 3, 4], [1, 2, 4]])
-      else
-        @test ray_indices(facets(S, Pos)) == IncidenceMatrix([[1, 3], [2, 3], [1, 2]])
-        @test vertex_and_ray_indices(facets(S, Pos)) ==
-          IncidenceMatrix([[1, 3, 4], [2, 3, 4], [1, 2, 4]])
-        @test IncidenceMatrix(facets(S, Pos)) ==
-          IncidenceMatrix([[1, 3, 4], [2, 3, 4], [1, 2, 4]])
-      end
-      @test vertex_indices(facets(S, Pos)) == IncidenceMatrix([[1], [1], [1]])
+      @test _check_im_perm_rows(ray_indices(facets(S, Pos)), [[2, 3], [1, 3], [1, 2]])
+      @test _check_im_perm_rows(vertex_and_ray_indices(facets(S, Pos)), [[2, 3, 4], [1, 3, 4], [1, 2, 4]])
+      @test _check_im_perm_rows(IncidenceMatrix(facets(S, Pos)), [[2, 3, 4], [1, 3, 4], [1, 2, 4]])
+      @test _check_im_perm_rows(vertex_indices(facets(S, Pos)), [[1], [1], [1]])
     end
-    @test facets(IncidenceMatrix, Pos) == IncidenceMatrix(
-      if T == QQFieldElem
+    @test _check_im_perm_rows(facets(IncidenceMatrix, Pos),
         [[2, 3, 4], [1, 3, 4], [1, 2, 4]]
-      else
-        [[1, 3, 4], [2, 3, 4], [1, 2, 4]]
-      end,
-    )
-    @test facet_indices(vertices(Pos)) == IncidenceMatrix([[1, 2, 3]])
-    @test IncidenceMatrix(vertices(Pos)) == IncidenceMatrix([[1, 2, 3]])
-    @test vertices(IncidenceMatrix, Pos) == IncidenceMatrix([[1, 2, 3]])
-    @test facet_indices(rays(Pos)) == (
-      if (T == QQFieldElem)
-        IncidenceMatrix([[2, 3], [1, 3], [1, 2]])
-      else
-        IncidenceMatrix([[1, 3], [2, 3], [1, 2]])
-      end
-    )
-    @test IncidenceMatrix(rays(Pos)) == (
-      if (T == QQFieldElem)
-        IncidenceMatrix([[2, 3], [1, 3], [1, 2]])
-      else
-        IncidenceMatrix([[1, 3], [2, 3], [1, 2]])
-      end
-    )
-    @test rays(IncidenceMatrix, Pos) == (
-      if (T == QQFieldElem)
-        IncidenceMatrix([[2, 3], [1, 3], [1, 2]])
-      else
-        IncidenceMatrix([[1, 3], [2, 3], [1, 2]])
-      end
     )
+    @test _check_im_perm_rows(facet_indices(vertices(Pos)), [[1, 2, 3]])
+    @test _check_im_perm_rows(IncidenceMatrix(vertices(Pos)), [[1, 2, 3]])
+    @test _check_im_perm_rows(vertices(IncidenceMatrix, Pos), [[1, 2, 3]])
+    @test _check_im_perm_rows(facet_indices(rays(Pos)), [[1, 3], [2, 3], [1, 2]])
+    @test _check_im_perm_rows(IncidenceMatrix(rays(Pos)), [[2, 3], [1, 3], [1, 2]])
+    @test _check_im_perm_rows(rays(IncidenceMatrix, Pos), [[1, 3], [2, 3], [1, 2]])
     @test facets(Pair, Pos) isa SubObjectIterator{Pair{Matrix{T},T}}
     @test facets(Pos) isa SubObjectIterator{AffineHalfspace{T}}
     @test facets(Halfspace, Pos) isa SubObjectIterator{AffineHalfspace{T}}
     @test affine_hull(point) isa SubObjectIterator{AffineHyperplane{T}}
-    @test affine_equation_matrix(affine_hull(point)) ==
-      matrix(f, [0 1 0 0; -1 0 1 0; 0 0 0 1])
-    @test Oscar.affine_matrix_for_polymake(affine_hull(point)) ==
-      [0 1 0 0; -1 0 1 0; 0 0 0 1]
+    @test affine_equation_matrix(affine_hull(point)) isa MatElem{T}
+    @test size(affine_equation_matrix(affine_hull(point))) == (3, 4)
+    @test Oscar.affine_matrix_for_polymake(affine_hull(point)) isa Polymake.Matrix{Oscar._scalar_type_to_polymake(T)}
+    @test size(Oscar.affine_matrix_for_polymake(affine_hull(point))) == (3, 4)
     @test length(affine_hull(point)) == 3
-    # TODO: restrict comparison to same scalar?
-    @test affine_hull(point) ==
-      [hyperplane(f, [1 0 0], 0), hyperplane(f, [0 1 0], 1), hyperplane(f, [0 0 1], 0)]
+    @test issetequal(affine_hull(point),
+                     [hyperplane(f, [1 0 0], 0), hyperplane(f, [0 1 0], 1), hyperplane(f, [0 0 1], 0)])
     @test n_facets(square) == 4
     @test lineality_dim(Q0) == 0
     @test n_rays(Q1) == 1
@@ -651,7 +617,7 @@
           @test halfspace_matrix_pair(facets(S, D)) == (A=matrix(R, vcat(A...)), b=b)
 
           @test ray_indices(facets(S, D)) == IncidenceMatrix(12, 0)
-          @test vertex_indices(facets(S, D)) == IncidenceMatrix([
+          @test _check_im_perm_rows(vertex_indices(facets(S, D)), [
             [1, 3, 5, 9, 10],
             [1, 2, 3, 4, 6],
             [1, 2, 5, 7, 8],
@@ -665,7 +631,7 @@
             [9, 10, 14, 17, 19],
             [15, 17, 18, 19, 20],
           ])
-          @test vertex_and_ray_indices(facets(S, D)) == IncidenceMatrix([
+          @test _check_im_perm_rows(vertex_and_ray_indices(facets(S, D)), [
             [1, 3, 5, 9, 10],
             [1, 2, 3, 4, 6],
             [1, 2, 5, 7, 8],
@@ -679,7 +645,7 @@
             [9, 10, 14, 17, 19],
             [15, 17, 18, 19, 20],
           ])
-          @test IncidenceMatrix(facets(S, D)) == IncidenceMatrix([
+          @test _check_im_perm_rows(IncidenceMatrix(facets(S, D)), [
             [1, 3, 5, 9, 10],
             [1, 2, 3, 4, 6],
             [1, 2, 5, 7, 8],
@@ -696,7 +662,7 @@
         end
       end
 
-      @test facets(IncidenceMatrix, D) == IncidenceMatrix([
+      @test _check_im_perm_rows(facets(IncidenceMatrix, D), [
         [1, 3, 5, 9, 10],
         [1, 2, 3, 4, 6],
         [1, 2, 5, 7, 8],
@@ -713,7 +679,7 @@
       @test facet_indices(rays(D)) == IncidenceMatrix(0, 12)
       @test IncidenceMatrix(rays(D)) == IncidenceMatrix(0, 12)
       @test rays(IncidenceMatrix, D) == IncidenceMatrix(0, 12)
-      @test facet_indices(vertices(D)) == IncidenceMatrix([
+      @test _check_im_perm_rows(facet_indices(vertices(D)), [
         [1, 2, 3],
         [2, 3, 6],
         [1, 2, 7],
@@ -735,7 +701,7 @@
         [9, 11, 12],
         [4, 9, 12],
       ])
-      @test IncidenceMatrix(vertices(D)) == IncidenceMatrix([
+      @test _check_im_perm_rows(IncidenceMatrix(vertices(D)), [
         [1, 2, 3],
         [2, 3, 6],
         [1, 2, 7],
@@ -757,7 +723,7 @@
         [9, 11, 12],
         [4, 9, 12],
       ])
-      @test vertices(IncidenceMatrix, D) == IncidenceMatrix([
+      @test _check_im_perm_rows(vertices(IncidenceMatrix, D), [
         [1, 2, 3],
         [2, 3, 6],
         [1, 2, 7],