diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl
index 898f936..380cee1 100644
--- a/engine-proto/Engine.Algebraic.jl
+++ b/engine-proto/Engine.Algebraic.jl
@@ -29,18 +29,6 @@ end
 dimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} =
   length(gens(base_ring(I))) - codimension(I, maxdepth)
 
-m2_ordering(R::MPolyRing) = Dict(
-  :lex => :Lex,
-  :deglex => :GLex,
-  :degrevlex => :GRevLex
-)[ordering(R)]
-
-string_m2(ring::MPolyRing) =
-  "QQ[$(join(symbols(ring), ", ")), MonomialOrder => $(m2_ordering(ring))]"
-
-string_m2(f::MPolyRingElem) =
-  replace(string(f), "//" => "/")
-
 # --- primitve elements ---
 
 abstract type Element{T} end
@@ -132,11 +120,11 @@ equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) -
 # --- constructions ---
 
 mutable struct Construction{T}
-  points::Vector{Point{T}}
-  spheres::Vector{Sphere{T}}
-  relations::Vector{Relation{T}}
+  points::Set{Point{T}}
+  spheres::Set{Sphere{T}}
+  relations::Set{Relation{T}}
   
-  function Construction{T}(; elements = Vector{Element{T}}(), relations = Vector{Relation{T}}()) where T
+  function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T
     allelements = union(elements, (rel.elements for rel in relations)...)
     new{T}(
       filter(elt -> isa(elt, Point), allelements),
@@ -161,10 +149,7 @@ function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T
   end
 end
 
-# output options:
-#  nothing - find a Gröbner basis
-#  :m2     - write a system of polynomials to a Macaulay2 file
-function realize(ctx::Construction{T}; output = nothing) where T
+function realize(ctx::Construction{T}) where T
   # collect coordinate names
   coordnamelist = Symbol[]
   eltenum = enumerate(Iterators.flatten((ctx.spheres, ctx.points)))
@@ -212,16 +197,8 @@ function realize(ctx::Construction{T}; output = nothing) where T
     push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts)
   end
   
-  if output == :m2
-    file = open("macaulay2/construction.m2", "w")
-    write(file, string(
-      "coordring = $(string_m2(coordring))\n",
-      "eqns = {\n  $(join(string_m2.(eqns), ",\n  "))\n}"
-    ))
-    close(file)
-  else
-    return (Generic.Ideal(coordring, eqns), eqns)
-  end
+  ## [test] (Generic.Ideal(coordring, eqns), eqns)
+  (nothing, eqns)
 end
 
 end
\ No newline at end of file
diff --git a/engine-proto/Engine.Numerical.jl b/engine-proto/Engine.Numerical.jl
index d1e14bd..48fb682 100644
--- a/engine-proto/Engine.Numerical.jl
+++ b/engine-proto/Engine.Numerical.jl
@@ -1,6 +1,5 @@
 module Numerical
 
-using Random: default_rng
 using LinearAlgebra
 using AbstractAlgebra
 using HomotopyContinuation:
@@ -29,16 +28,16 @@ end
 
 # --- sampling ---
 
-function real_samples(F::AbstractSystem, dim; rng = default_rng())
+function real_samples(F::AbstractSystem, dim)
   # 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(rng, nvariables(F))) for _ in 1:dim)...
+    (normalize(randn(nvariables(F))) for _ in 1:dim)...
   ))
   cut = LinearSubspace(normals, fill(0., dim))
-  filter(isreal, results(witness_set(F, cut, seed = 0x1974abba)))
+  filter(isreal, results(witness_set(F, cut)))
 end
 
 AbstractAlgebra.evaluate(pt::Point, vals::Vector{<:RingElement}) =
diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl
index af4128d..49011c6 100644
--- a/engine-proto/Engine.jl
+++ b/engine-proto/Engine.jl
@@ -23,55 +23,106 @@ using GLMakie
 
 CoeffType = Rational{Int64}
 
-spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
-tangencies = [
-  Engine.AlignsWithBy{CoeffType}(
-    spheres[n],
-    spheres[mod1(n+1, length(spheres))],
-    CoeffType(1)
-  )
-  for n in 1:3
-]
-ctx_tan_sph = Engine.Construction{CoeffType}(elements = spheres, relations = tangencies)
+##a = Engine.Point{CoeffType}()
+##s = Engine.Sphere{CoeffType}()
+##a_on_s = Engine.LiesOn{CoeffType}(a, s)
+##ctx = Engine.Construction{CoeffType}(elements = Set([a]), relations= Set([a_on_s]))
+##ideal_a_s = Engine.realize(ctx)
+##println("A point on a sphere: $(Engine.dimension(ideal_a_s)) degrees of freedom")
+
+##b = Engine.Point{CoeffType}()
+##b_on_s = Engine.LiesOn{CoeffType}(b, s)
+##Engine.push!(ctx, b)
+##Engine.push!(ctx, s)
+##Engine.push!(ctx, b_on_s)
+##ideal_ab_s, eqns_ab_s = Engine.realize(ctx)
+##freedom = Engine.dimension(ideal_ab_s)
+##println("Two points on a sphere: $freedom degrees of freedom")
+
+##spheres = [Engine.Sphere{CoeffType}() for _ in 1:3]
+##tangencies = [
+##  Engine.AlignsWithBy{CoeffType}(
+##    spheres[n],
+##    spheres[mod1(n+1, length(spheres))],
+##    CoeffType(-1//1)
+##  )
+##  for n in 1:3
+##]
+##tangencies = [
+  ##Engine.LiesOn{CoeffType}(points[1], spheres[2]),
+  ##Engine.LiesOn{CoeffType}(points[1], spheres[3]),
+  ##Engine.LiesOn{CoeffType}(points[2], spheres[3]),
+  ##Engine.LiesOn{CoeffType}(points[2], spheres[1]),
+  ##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))
 ##ideal_tan_sph, eqns_tan_sph = Engine.realize(ctx_tan_sph)
-Engine.realize(ctx_tan_sph, output = :m2)
 ##freedom = Engine.dimension(ideal_tan_sph)
 ##println("Three mutually tangent spheres: $freedom degrees of freedom")
 
+points = [Engine.Point{CoeffType}() for _ in 1:3]
+spheres = [Engine.Sphere{CoeffType}() for _ in 1:2]
+ctx_joined = Engine.Construction{CoeffType}(
+  elements = Set([points; spheres]),
+  relations= Set([
+    Engine.LiesOn{CoeffType}(pt, sph)
+    for pt in points for sph in spheres
+  ])
+)
+ideal_joined, eqns_joined = Engine.realize(ctx_joined)
+freedom = Engine.dimension(ideal_joined)
+println("$(length(points)) points on $(length(spheres)) spheres: $freedom degrees of freedom")
+
 # --- test rational cut ---
 
-##coordring = base_ring(ideal_tan_sph)
-##vbls = Variable.(symbols(coordring))
+coordring = base_ring(ideal_joined)
+vbls = Variable.(symbols(coordring))
 
 # test a random witness set
-##system = CompiledSystem(System(eqns_tan_sph, variables = vbls))
-##norm2 = vec -> real(dot(conj.(vec), vec))
-##rng = MersenneTwister(6071)
-##n_planes = 6
-##samples = []
-##for _ in 1:n_planes
-##  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)
-##    end
-##  end
-##end
-##println("Found $(length(samples)) sample solutions")
+system = CompiledSystem(System(eqns_joined, variables = vbls))
+norm2 = vec -> real(dot(conj.(vec), vec))
+Random.seed!(6071)
+n_planes = 3
+samples = []
+for _ in 1:n_planes
+  real_solns = solution.(Engine.Numerical.real_samples(system, freedom))
+  for soln in real_solns
+    if all(norm2(soln - samp) > 1e-4*length(gens(coordring)) for samp in samples)
+      push!(samples, soln)
+    end
+  end
+end
+println("$(length(samples)) sample solutions:")
+for soln in samples
+  ## display([vbls round.(soln, digits = 6)]) ## [verbose]
+  k_sq = abs2(soln[1])
+  if abs2(soln[end-2]) > 1e-12
+    if k_sq < 1e-12
+      println("  center at infinity: z coordinates $(round(soln[end], digits = 6)) and $(round(soln[end-1], digits = 6))")
+    else
+      sum_sq = soln[4]^2 + soln[7]^2 + soln[end-2]^2 / k_sq
+      println("  center on z axis:   r² = $(round(1/k_sq, digits = 6)), x² + y² + h² = $(round(sum_sq, digits = 6))")
+    end
+  else
+    sum_sq = sum(soln[[4, 7, 10]] .^ 2)
+    println("  center at origin:   r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))")
+  end
+end
 
 # show a sample solution
-##function show_solution(ctx, vals)
-##  # evaluate elements
-##  real_vals = real.(vals)
-##  disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points]
-##  disp_spheres = [Engine.Numerical.evaluate(sph, real_vals) for sph in ctx.spheres]
-##  
-##  # create scene
-##  scene = Scene()
-##  cam3d!(scene)
-##  scatter!(scene, disp_points, color = :green)
-##  for sph in disp_spheres
-##    mesh!(scene, sph, color = :gray)
-##  end
-##  scene
-##end
\ No newline at end of file
+function show_solution(ctx, vals)
+  # evaluate elements
+  real_vals = real.(vals)
+  disp_points = [Engine.Numerical.evaluate(pt, real_vals) for pt in ctx.points]
+  disp_spheres = [Engine.Numerical.evaluate(sph, real_vals) for sph in ctx.spheres]
+  
+  # create scene
+  scene = Scene()
+  cam3d!(scene)
+  scatter!(scene, disp_points, color = :green)
+  for sph in disp_spheres
+    mesh!(scene, sph, color = :gray)
+  end
+  scene
+end
\ No newline at end of file
diff --git a/engine-proto/macaulay2/engine.m2 b/engine-proto/macaulay2/engine.m2
deleted file mode 100644
index 43eb407..0000000
--- a/engine-proto/macaulay2/engine.m2
+++ /dev/null
@@ -1,3 +0,0 @@
-needsPackage "TriangularSets"
-
-mprod = (v, w) -> (v#0*w#1 + w#0*v#1) / 2 - v#2*w#2 - v#3*w#3 - v#4*w#4
\ No newline at end of file