From 291d5c8ff685b642fd3535b6c46179c9f134c469 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Thu, 15 Feb 2024 13:28:01 -0800 Subject: [PATCH] Study mutually tangent spheres with two fixed --- engine-proto/Engine.Algebraic.jl | 24 ++++----- engine-proto/Engine.jl | 93 +++++++++++++++++--------------- 2 files changed, 62 insertions(+), 55 deletions(-) diff --git a/engine-proto/Engine.Algebraic.jl b/engine-proto/Engine.Algebraic.jl index 380cee1..b6fc7a7 100644 --- a/engine-proto/Engine.Algebraic.jl +++ b/engine-proto/Engine.Algebraic.jl @@ -186,19 +186,19 @@ 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 - ## [test] (Generic.Ideal(coordring, eqns), eqns) - (nothing, eqns) + (Generic.Ideal(coordring, eqns), eqns) + ## [test] (nothing, eqns) end end \ No newline at end of file diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 49011c6..0410f6d 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -39,15 +39,15 @@ CoeffType = Rational{Int64} ##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 -##] +spheres = [Engine.Sphere{CoeffType}() for _ in 1:3] +tangencies = [ + Engine.AlignsWithBy{CoeffType}( + spheres[n], + spheres[mod1(n+1, length(spheres))], + CoeffType(-1)^n + ) + for n in 1:3 +] ##tangencies = [ ##Engine.LiesOn{CoeffType}(points[1], spheres[2]), ##Engine.LiesOn{CoeffType}(points[1], spheres[3]), @@ -56,34 +56,41 @@ CoeffType = Rational{Int64} ##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) -##freedom = Engine.dimension(ideal_tan_sph) -##println("Three mutually tangent spheres: $freedom degrees of freedom") +ctx_tan_sph = Engine.Construction{CoeffType}(elements = Set(spheres), relations = Set(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) +println("Three mutually tangent spheres, with two fixed: $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") +##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_joined) +coordring = base_ring(small_ideal_tan_sph) vbls = Variable.(symbols(coordring)) # test a random witness set -system = CompiledSystem(System(eqns_joined, variables = vbls)) +system = CompiledSystem(System(small_eqns_tan_sph, variables = vbls)) norm2 = vec -> real(dot(conj.(vec), vec)) Random.seed!(6071) -n_planes = 3 +n_planes = 36 samples = [] for _ in 1:n_planes real_solns = solution.(Engine.Numerical.real_samples(system, freedom)) @@ -94,21 +101,21 @@ for _ in 1:n_planes 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 +##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)