engine-proto/Engine.jl

@@ -263,67 +263,54 @@ println("Two points on a sphere: ", freedom, " degrees of freedom")
 
 # --- test rational cut ---
 
-cut_coeffs = [
+coordring = base_ring(ideal_ab_s)
-  1 1 1 0 0 0 1 1 1 1 1;
+vbls = Variable.(symbols(coordring))
-  2 1 1 0 0 0 1 2 1 1 1;
+##cut_system = CompiledSystem(System([eqns_ab_s; cut], variables = vbls))
-  1 2 0 0 0 0 1 1 0 1 2
+##cut_result = HomotopyContinuation.solve(cut_system)
-]
+##println("non-singular solutions:")
-cut = [
+##for soln in solutions(cut_result)
-  sum(vcat(cf[1:3] .* a.coords, cf[4:6] .* b.coords, cf[7:end] .* (s.coords - [0, 0, 0, 0, 1])))
+##  display(soln)
-  for cf in eachrow(cut_coeffs)
+##end
-]
+##println("singular solutions:")
-cut_ideal_ab_s = Generic.Ideal(base_ring(ideal_ab_s), [gens(ideal_ab_s); cut])
+##for sing in singular(cut_result)
-cut_freedom = Engine.dimension(cut_ideal_ab_s)
+##  display(sing.solution)
-println("Two points on a sphere, after cut: ", cut_freedom, " degrees of freedom")
+##end
-if cut_freedom == 0
+
-  coordring = base_ring(ideal_ab_s)
+# test a random witness set
-  vbls = Variable.(symbols(coordring))
+system = CompiledSystem(System(eqns_ab_s, variables = vbls))
-  cut_system = System([eqns_ab_s; cut], variables = vbls)
+max_slope = 2
-  ##cut_result = HomotopyContinuation.solve(cut_system)
+binom = Binomial(2max_slope, 1/2)
-  ##println("non-singular solutions:")
+Random.seed!(6071)
-  ##for soln in solutions(cut_result)
+samples = []
-  ##  display(soln)
+for _ in 1:3
-  ##end
+  cut_matrix = rand(binom, freedom, length(gens(coordring))) .- max_slope
-  ##println("singular solutions:")
+  ##cut_matrix = [
-  ##for sing in singular(cut_result)
+  ##  1 1 1 1 0 1 1 0 1 1 0;
-  ##  display(sing.solution)
+  ##  1 2 1 2 0 1 1 0 1 1 0;
-  ##end
+  ##  1 1 0 1 0 1 2 0 2 0 0
-  
+  ##]
-  # test a random witness set
+  sph_z_ind = indexin([sph.coords[5] for sph in ctx.spheres], gens(coordring))
-  max_slope = 2
+  cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)]
-  binom = Binomial(2max_slope, 1/2)
+  println("sphere z variables: ", vbls[sph_z_ind])
-  Random.seed!(6071)
+  display(cut_matrix)
-  samples = []
+  display(cut_offset)
-  for _ in 1:3
+  cut_subspace = LinearSubspace(cut_matrix, cut_offset)
-    cut_matrix = rand(binom, freedom, length(gens(coordring))) .- max_slope
+  wtns = witness_set(system, cut_subspace)
-    ##cut_matrix = [
+  append!(samples, solution.(filter(isreal, results(wtns))))
-    ##  1 1 1 1 0 1 1 0 1 1 0;
+end
-    ##  1 2 1 2 0 1 1 0 1 1 0;
+println("witness solutions:")
-    ##  1 1 0 1 0 1 2 0 2 0 0
+for soln in samples
-    ##]
+  display([vbls round.(soln, digits = 6)])
-    sph_z_ind = indexin([sph.coords[5] for sph in ctx.spheres], gens(coordring))
+  k_sq = abs2(soln[1])
-    cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)]
+  if abs2(soln[end-2]) > 1e-12
-    println("sphere z variables: ", vbls[sph_z_ind])
+    if k_sq < 1e-12
-    display(cut_matrix)
+      println("center at infinity: z coordinates $(round(soln[end], digits = 6)) and $(round(soln[end-1], digits = 6))}")
-    display(cut_offset)
-    cut_subspace = LinearSubspace(cut_matrix, cut_offset)
-    wtns = witness_set(System(eqns_ab_s, variables = vbls), cut_subspace)
-    append!(samples, solution.(filter(isreal, results(wtns))))
-  end
-  println("witness solutions:")
-  for soln in samples
-    display([vbls round.(soln, digits = 6)])
-    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)
+      sum_sq = soln[4]^2 + soln[7]^2 + soln[end-2]^2 / k_sq
-      println("center at origin: r² = $(round(1/k_sq, digits = 6)); x² + y² + z² = $(round(sum_sq, digits = 6))")
+      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