From ba365174d3a643c033a04178706fe04bfe2bc555 Mon Sep 17 00:00:00 2001
From: Aaron Fenyes <aaron.fenyes@fareycircles.ooo>
Date: Thu, 15 Feb 2024 16:16:06 -0800
Subject: [PATCH] Find real solutions for three mutually tangent spheres

I'm not sure why the solver wasn't working before. It might've been just
an unlucky random number draw.
---
 engine-proto/Engine.Algebraic.jl | 20 ++++++++++----------
 engine-proto/Engine.jl           | 18 +++++++++---------
 2 files changed, 19 insertions(+), 19 deletions(-)

diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl
index f14f58c..2d28017 100644
--- a/engine-proto/Engine.Algebraic.jl
+++ b/engine-proto/Engine.Algebraic.jl
@@ -186,16 +186,16 @@ function realize(ctx::Construction{T}) where T
   # add relations to center, orient, and scale the construction
   # [to do] the scaling constraint, as written, can be impossible to satisfy
   #         when all of the spheres have to go through the origin
-  ##if !isempty(ctx.points)
-  ##  append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3])
-  ##end
-  ##if !isempty(ctx.spheres)
-  ##  append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4])
-  ##end
-  ##n_elts = length(ctx.points) + length(ctx.spheres)
-  ##if n_elts > 0
-  ##  push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts)
-  ##end
+  if !isempty(ctx.points)
+    append!(eqns, [sum(pt.coords[k] for pt in ctx.points) for k in 1:3])
+  end
+  if !isempty(ctx.spheres)
+    append!(eqns, [sum(sph.coords[k] for sph in ctx.spheres) for k in 3:4])
+  end
+  n_elts = length(ctx.points) + length(ctx.spheres)
+  if n_elts > 0
+    push!(eqns, sum(elt.vec[2] for elt in Iterators.flatten((ctx.points, ctx.spheres))) - n_elts)
+  end
   
   (Generic.Ideal(coordring, eqns), eqns)
   ## [test] (nothing, eqns)
diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl
index 7146e80..e92eed4 100644
--- a/engine-proto/Engine.jl
+++ b/engine-proto/Engine.jl
@@ -59,13 +59,13 @@ 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,11 +83,11 @@ 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))
 rng = MersenneTwister(6071)
 n_planes = 3