diff --git a/cython/pycddlib.pxi b/cython/pycddlib.pxi
index 4f1c0b7..1dae9e9 100644
--- a/cython/pycddlib.pxi
+++ b/cython/pycddlib.pxi
@@ -78,9 +78,9 @@ cdef _tmpread(libc.stdio.FILE *pfile):
 cdef _get_set(set_type set_):
     # create Python Set from given set_type
     cdef unsigned long elem
-    return frozenset(
+    return {
         elem for elem from 0 <= elem < set_[0] if set_member(elem + 1, set_)
-    )
+    }
 
 cdef _set_set(set_type set_, pset):
     # set elements of set_type by elements from a Python Container
@@ -99,11 +99,11 @@ cdef _get_dd_setfam(dd_SetFamilyPtr setfam):
     if setfam == NULL:
         raise ValueError("failed to get set family")
     result = [
-        frozenset(
+        {
             elem
             for elem from 0 <= elem < setfam.setsize
             if set_member(elem + 1, setfam.set[i])
-        )
+        }
         for i from 0 <= i < setfam.famsize
     ]
     dd_FreeSetFamily(setfam)
diff --git a/docs/source/matrix.rst b/docs/source/matrix.rst
index 82e19d8..fa4aa83 100644
--- a/docs/source/matrix.rst
+++ b/docs/source/matrix.rst
@@ -2,6 +2,7 @@
 
     import cdd.gmp
     from fractions import Fraction
+    from pprint import pprint
 
 Sets of Linear Inequalities and Generators
 ==========================================
@@ -9,40 +10,20 @@ Sets of Linear Inequalities and Generators
 Declaring matrices, and checking some attributes:
 
 >>> mat1 = cdd.gmp.matrix_from_array([[1, 2],[3, 4]])
->>> print(mat1) # doctest: +NORMALIZE_WHITESPACE
-begin
- 2 2 rational
- 1 2
- 3 4
-end
->>> print(mat1.array)
+>>> pprint(mat1.array)
 [[Fraction(1, 1), Fraction(2, 1)], [Fraction(3, 1), Fraction(4, 1)]]
 >>> cdd.gmp.matrix_append_to(mat1, cdd.gmp.matrix_from_array([[5,6]]))
->>> print(mat1) # doctest: +NORMALIZE_WHITESPACE
-begin
- 3 2 rational
- 1 2
- 3 4
- 5 6
-end
->>> print(mat1.array)
-[[Fraction(1, 1), Fraction(2, 1)], [Fraction(3, 1), Fraction(4, 1)], [Fraction(5, 1), Fraction(6, 1)]]
+>>> pprint(mat1.array)
+[[Fraction(1, 1), Fraction(2, 1)],
+ [Fraction(3, 1), Fraction(4, 1)],
+ [Fraction(5, 1), Fraction(6, 1)]]
 
 Canonicalizing:
 
->>> mat = cdd.gmp.matrix_from_array([[2, 1, 2, 3], [0, 1, 2, 3], [3, 0, 1, 2], [0, -2, -4, -6]])
->>> cdd.gmp.matrix_canonicalize(mat)  # oops... must specify rep_type!
-Traceback (most recent call last):
-    ...
-ValueError: rep_type unspecified
->>> mat.rep_type = cdd.RepType.INEQUALITY
+>>> array = [[2, 1, 2, 3], [0, 1, 2, 3], [3, 0, 1, 2], [0, -2, -4, -6]]
+>>> mat = cdd.gmp.matrix_from_array(array, rep_type=cdd.RepType.INEQUALITY)
 >>> cdd.gmp.matrix_canonicalize(mat)
-(frozenset({1, 3}), frozenset({0}))
->>> print(mat) # doctest: +NORMALIZE_WHITESPACE
-H-representation
-linearity 1  1
-begin
- 2 4 rational
- 0 1 2 3
- 3 0 1 2
-end
+({1, 3}, {0})
+>>> pprint(mat.array)
+[[Fraction(0, 1), Fraction(1, 1), Fraction(2, 1), Fraction(3, 1)],
+ [Fraction(3, 1), Fraction(0, 1), Fraction(1, 1), Fraction(2, 1)]]
diff --git a/docs/source/polyhedron.rst b/docs/source/polyhedron.rst
index 19ffb89..3dc2589 100644
--- a/docs/source/polyhedron.rst
+++ b/docs/source/polyhedron.rst
@@ -1,35 +1,26 @@
 .. testsetup::
 
-   import cdd
+    import cdd
+    from pprint import pprint
 
 Working With Polyhedron Representations
 =======================================
 
 This is the sampleh1.ine example that comes with cddlib.
 
->>> mat = cdd.matrix_from_array([[2, -1, -1, 0],[0, 1, 0, 0],[0, 0, 1, 0]])
->>> mat.rep_type = cdd.RepType.INEQUALITY
+>>> array = [[2, -1, -1, 0], [0, 1, 0, 0], [0, 0, 1, 0]]
+>>> mat = cdd.matrix_from_array(array, rep_type=cdd.RepType.INEQUALITY)
 >>> poly = cdd.polyhedron_from_matrix(mat)
->>> print(poly) # doctest: +NORMALIZE_WHITESPACE
-begin
- 3 4 real
- 2 -1 -1 0
- 0 1 0 0
- 0 0 1 0
-end
 >>> ext = cdd.copy_generators(poly)
->>> print(ext) # doctest: +NORMALIZE_WHITESPACE
-V-representation
-linearity 1  4
-begin
- 4 4 real
- 1 0 0 0
- 1 2 0 0
- 1 0 2 0
- 0 0 0 1
-end
->>> print(list(ext.lin_set)) # note: first row is 0, so fourth row is 3
-[3]
+>>> ext.rep_type
+<RepType.GENERATOR: 2>
+>>> pprint(ext.array) # doctest: +NORMALIZE_WHITESPACE
+[[1.0, 0.0, 0.0, 0.0],
+ [1.0, 2.0, 0.0, 0.0],
+ [1.0, 0.0, 2.0, 0.0],
+ [0.0, 0.0, 0.0, 1.0]]
+>>> ext.lin_set # note: first row is 0, so fourth row is 3
+{3}
 
 
 The following example illustrates how to get adjacencies and incidences.
@@ -44,15 +35,12 @@ The following example illustrates how to get adjacencies and incidences.
 >>> poly = cdd.polyhedron_from_matrix(mat)
 >>> # The V-representation can be printed in the usual way:
 >>> gen = cdd.copy_generators(poly)
->>> print(gen) # doctest: +NORMALIZE_WHITESPACE
-V-representation
-begin
- 4 3 real
- 1 1 -1
- 1 1 1
- 1 -1 1
- 1 -1 -1
-end
+>>> gen.rep_type
+<RepType.GENERATOR: 2>
+>>> pprint(gen.array)
+[[1.0, 1.0, -1.0], [1.0, 1.0, 1.0], [1.0, -1.0, 1.0], [1.0, -1.0, -1.0]]
+>>> gen.lin_set
+set()
 >>> # graphical depiction of vertices and faces:
 >>> #
 >>> #   2---(3)---1
@@ -67,39 +55,40 @@ end
 >>> # vertex 1 is adjacent to vertices 0 and 2
 >>> # vertex 2 is adjacent to vertices 1 and 3
 >>> # vertex 3 is adjacent to vertices 0 and 2
->>> print([list(x) for x in cdd.copy_adjacency(poly)])
-[[1, 3], [0, 2], [1, 3], [0, 2]]
+>>> cdd.copy_adjacency(poly)
+[{1, 3}, {0, 2}, {1, 3}, {0, 2}]
 >>> # vertex 0 is the intersection of faces (1) and (2)
 >>> # vertex 1 is the intersection of faces (2) and (3)
 >>> # vertex 2 is the intersection of faces (0) and (3)
 >>> # vertex 3 is the intersection of faces (0) and (1)
->>> print([list(x) for x in cdd.copy_incidence(poly)])
-[[1, 2], [2, 3], [0, 3], [0, 1]]
+>>> cdd.copy_incidence(poly)
+[{1, 2}, {2, 3}, {0, 3}, {0, 1}]
 >>> # face (0) is adjacent to faces (1) and (3)
 >>> # face (1) is adjacent to faces (0) and (2)
 >>> # face (2) is adjacent to faces (1) and (3)
 >>> # face (3) is adjacent to faces (0) and (2)
->>> print([list(x) for x in cdd.copy_input_adjacency(poly)])
-[[1, 3], [0, 2], [1, 3], [0, 2], []]
+>>> cdd.copy_input_adjacency(poly)
+[{1, 3}, {0, 2}, {1, 3}, {0, 2}, set()]
 >>> # face (0) intersects with vertices 2 and 3
 >>> # face (1) intersects with vertices 0 and 3
 >>> # face (2) intersects with vertices 0 and 1
 >>> # face (3) intersects with vertices 1 and 2
->>> print([list(x) for x in cdd.copy_input_incidence(poly)])
-[[2, 3], [0, 3], [0, 1], [1, 2], []]
+>>> cdd.copy_input_incidence(poly)
+[{2, 3}, {0, 3}, {0, 1}, {1, 2}, set()]
 >>> # add a vertex, and construct new polyhedron
 >>> cdd.matrix_append_to(gen, cdd.matrix_from_array([[1, 0, 2]]))
 >>> vpoly = cdd.polyhedron_from_matrix(gen)
->>> print(cdd.copy_inequalities(vpoly)) # doctest: +NORMALIZE_WHITESPACE
-H-representation
-begin
- 5 3 real
- 1 0 1
- 2 1 -1
- 1 1 0
- 2 -1 -1
- 1 -1 0
-end
+>>> vmat = cdd.copy_inequalities(vpoly)
+>>> vmat.rep_type
+<RepType.INEQUALITY: 1>
+>>> pprint(vmat.array)
+[[1.0, 0.0, 1.0],
+ [2.0, 1.0, -1.0],
+ [1.0, 1.0, 0.0],
+ [2.0, -1.0, -1.0],
+ [1.0, -1.0, 0.0]]
+>>> vmat.lin_set
+set()
 >>> # so now we have:
 >>> # 0 <= 1 + x2
 >>> # 0 <= 2 + x1 - x2
@@ -123,14 +112,14 @@ end
 >>> #   3---(0)---0
 >>> #
 >>> # for each face, list adjacent faces
->>> print([list(x) for x in cdd.copy_adjacency(vpoly)])
-[[2, 4], [2, 3], [0, 1], [1, 4], [0, 3]]
+>>> cdd.copy_adjacency(vpoly)
+[{2, 4}, {2, 3}, {0, 1}, {1, 4}, {0, 3}]
 >>> # for each face, list adjacent vertices
->>> print([list(x) for x in cdd.copy_incidence(vpoly)])
-[[0, 3], [2, 4], [2, 3], [1, 4], [0, 1]]
+>>> cdd.copy_incidence(vpoly)
+[{0, 3}, {2, 4}, {2, 3}, {1, 4}, {0, 1}]
 >>> # for each vertex, list adjacent vertices
->>> print([list(x) for x in cdd.copy_input_adjacency(vpoly)])
-[[1, 3], [0, 4], [3, 4], [0, 2], [1, 2]]
+>>> cdd.copy_input_adjacency(vpoly)
+[{1, 3}, {0, 4}, {3, 4}, {0, 2}, {1, 2}]
 >>> # for each vertex, list adjacent faces
->>> print([list(x) for x in cdd.copy_input_incidence(vpoly)])
-[[0, 4], [3, 4], [1, 2], [0, 2], [1, 3]]
+>>> cdd.copy_input_incidence(vpoly)
+[{0, 4}, {3, 4}, {1, 2}, {0, 2}, {1, 3}]
diff --git a/test/test_adjacency_list.py b/test/test_adjacency_list.py
index bde720d..ae4a7eb 100644
--- a/test/test_adjacency_list.py
+++ b/test/test_adjacency_list.py
@@ -32,16 +32,16 @@ def test_make_vertex_adjacency_list() -> None:
     # The first vertex is adjacent to the second, fourth and eighth
     # (note the conversion to a pythonic numbering system)
     adjacency_list = [
-        [1, 3, 7],
-        [0, 2, 6],
-        [1, 3, 4],
-        [0, 2, 5],
-        [2, 5, 6],
-        [3, 4, 7],
-        [1, 4, 7],
-        [0, 5, 6],
+        {1, 3, 7},
+        {0, 2, 6},
+        {1, 3, 4},
+        {0, 2, 5},
+        {2, 5, 6},
+        {3, 4, 7},
+        {1, 4, 7},
+        {0, 5, 6},
     ]
-    assert adjacency == [frozenset(x) for x in adjacency_list]
+    assert adjacency == adjacency_list
 
 
 def test_make_facet_adjacency_list() -> None:
@@ -60,14 +60,14 @@ def test_make_facet_adjacency_list() -> None:
     poly = cdd.polyhedron_from_matrix(mat)
 
     adjacency_list = [
-        [1, 2, 3, 4, 6],
-        [0, 2, 3, 5],
-        [0, 1, 4, 5],
-        [0, 1, 5, 6],
-        [0, 2, 5, 6],
-        [1, 2, 3, 4, 6],
-        [0, 3, 4, 5],
+        {1, 2, 3, 4, 6},
+        {0, 2, 3, 5},
+        {0, 1, 4, 5},
+        {0, 1, 5, 6},
+        {0, 2, 5, 6},
+        {1, 2, 3, 4, 6},
+        {0, 3, 4, 5},
     ]
 
     adjacency = cdd.copy_input_adjacency(poly)
-    assert adjacency == [frozenset(x) for x in adjacency_list]
+    assert adjacency == adjacency_list
diff --git a/test/test_incidence.py b/test/test_incidence.py
index 35959e5..e02da7a 100644
--- a/test/test_incidence.py
+++ b/test/test_incidence.py
@@ -39,7 +39,7 @@ def test_vertex_incidence_cube() -> None:
         {0, 1, 5},
         {0, 1, 2},
     ]
-    assert incidence == [frozenset(x) for x in incidence_list]
+    assert incidence == incidence_list
 
 
 def test_vertex_incidence_vtest_vo() -> None:
@@ -72,7 +72,7 @@ def test_vertex_incidence_vtest_vo() -> None:
     ]
 
     incidence = cdd.copy_incidence(poly)
-    assert incidence == [frozenset(x) for x in incidence_list]
+    assert incidence == incidence_list
 
 
 def test_facet_incidence_cube() -> None:
@@ -111,7 +111,7 @@ def test_facet_incidence_cube() -> None:
         {1, 2, 4, 6},
         set(),
     ]
-    assert incidence == [frozenset(x) for x in incidence_list]
+    assert incidence == incidence_list
 
 
 def test_facet_incidence_vtest_vo() -> None:
@@ -140,4 +140,4 @@ def test_facet_incidence_vtest_vo() -> None:
         {0, 4, 7, 8},
     ]
 
-    assert cdd.copy_input_incidence(poly) == [frozenset(x) for x in incidence_list]
+    assert cdd.copy_input_incidence(poly) == incidence_list