diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl
index b6fc7a7..f14f58c 100644
--- a/engine-proto/Engine.Algebraic.jl
+++ b/engine-proto/Engine.Algebraic.jl
@@ -120,11 +120,11 @@ equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) -
 # --- constructions ---
 
 mutable struct Construction{T}
-  points::Set{Point{T}}
-  spheres::Set{Sphere{T}}
-  relations::Set{Relation{T}}
+  points::Vector{Point{T}}
+  spheres::Vector{Sphere{T}}
+  relations::Vector{Relation{T}}
   
-  function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T
+  function Construction{T}(; elements = Vector{Element{T}}(), relations = Vector{Relation{T}}()) where T
     allelements = union(elements, (rel.elements for rel in relations)...)
     new{T}(
       filter(elt -> isa(elt, Point), allelements),
diff --git a/engine-proto/Engine.Numerical.jl b/engine-proto/Engine.Numerical.jl
index 48fb682..c54e71c 100644
--- a/engine-proto/Engine.Numerical.jl
+++ b/engine-proto/Engine.Numerical.jl
@@ -1,5 +1,6 @@
 module Numerical
 
+using Random: default_rng
 using LinearAlgebra
 using AbstractAlgebra
 using HomotopyContinuation:
@@ -28,16 +29,18 @@ end
 
 # --- sampling ---
 
-function real_samples(F::AbstractSystem, dim)
+function real_samples(F::AbstractSystem, dim; rng = default_rng())
   # choose a random real hyperplane of codimension `dim` by intersecting
   # hyperplanes whose normal vectors are uniformly distributed over the unit
   # sphere
   # [to do] guard against the unlikely event that one of the normals is zero
-  normals = transpose(hcat(
-    (normalize(randn(nvariables(F))) for _ in 1:dim)...
-  ))
-  cut = LinearSubspace(normals, fill(0., dim))
-  filter(isreal, results(witness_set(F, cut)))
+  ##normals = transpose(hcat(
+  ##  (normalize(randn(rng, nvariables(F))) for _ in 1:dim)...
+  ##))
+  ##cut = LinearSubspace(normals, fill(0., dim))
+  ##filter(isreal, results(witness_set(F, cut, seed = 0x8af341df)))
+  ##filter(isreal, results(witness_set(F, seed = 0x8af341df)))
+  results(witness_set(F, seed = 0x8af341df))
 end
 
 AbstractAlgebra.evaluate(pt::Point, vals::Vector{<:RingElement}) =
diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl
index 0410f6d..9f7d150 100644
--- a/engine-proto/Engine.jl
+++ b/engine-proto/Engine.jl
@@ -56,16 +56,16 @@ tangencies = [
   ##Engine.LiesOn{CoeffType}(points[3], spheres[1]),
   ##Engine.LiesOn{CoeffType}(points[3], spheres[2])
 ##]
-ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(tangencies))
+ctx_tan_sph = Engine.Construction{CoeffType}(elements = spheres, relations = tangencies)
 ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph)
 ##small_eqns_tan_sph = eqns_tan_sph
-small_eqns_tan_sph = [
-  eqns_tan_sph;
-  spheres[2].coords - [1, 0, 0, 0, 1];
-  spheres[3].coords - [1, 0, 0, 0, -1];
-]
-small_ideal_tan_sph = Generic.Ideal(base_ring(ideal_tan_sph), small_eqns_tan_sph)
-freedom = Engine.dimension(small_ideal_tan_sph)
+##small_eqns_tan_sph = [
+##  eqns_tan_sph;
+##  spheres[2].coords - [1, 0, 0, 0, 1];
+##  spheres[3].coords - [1, 0, 0, 0, -1];
+##]
+##small_ideal_tan_sph = Generic.Ideal(base_ring(ideal_tan_sph), small_eqns_tan_sph)
+freedom = Engine.dimension(ideal_tan_sph)
 println("Three mutually tangent spheres, with two fixed: $freedom degrees of freedom")
 
 ##points = [Engine.Point{CoeffType}() for _ in 1:3]
@@ -83,17 +83,17 @@ println("Three mutually tangent spheres, with two fixed: $freedom degrees of fre
 
 # --- test rational cut ---
 
-coordring = base_ring(small_ideal_tan_sph)
+coordring = base_ring(ideal_tan_sph)
 vbls = Variable.(symbols(coordring))
 
 # test a random witness set
-system = CompiledSystem(System(small_eqns_tan_sph, variables = vbls))
+system = CompiledSystem(System(eqns_tan_sph, variables = vbls))
 norm2 = vec -> real(dot(conj.(vec), vec))
-Random.seed!(6071)
-n_planes = 36
+rng = MersenneTwister(6701)
+n_planes = 6
 samples = []
 for _ in 1:n_planes
-  real_solns = solution.(Engine.Numerical.real_samples(system, freedom))
+  real_solns = solution.(Engine.Numerical.real_samples(system, freedom, rng = rng))
   for soln in real_solns
     if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples)
       push!(samples, soln)