diff --git a/engine-proto/Engine.jl b/engine-proto/Engine.jl index 546bf21..4bce0d7 100644 --- a/engine-proto/Engine.jl +++ b/engine-proto/Engine.jl @@ -263,67 +263,54 @@ println("Two points on a sphere: ", freedom, " degrees of freedom") # --- test rational cut --- -cut_coeffs = [ - 1 1 1 0 0 0 1 1 1 1 1; - 2 1 1 0 0 0 1 2 1 1 1; - 1 2 0 0 0 0 1 1 0 1 2 -] -cut = [ - sum(vcat(cf[1:3] .* a.coords, cf[4:6] .* b.coords, cf[7:end] .* (s.coords - [0, 0, 0, 0, 1]))) - for cf in eachrow(cut_coeffs) -] -cut_ideal_ab_s = Generic.Ideal(base_ring(ideal_ab_s), [gens(ideal_ab_s); cut]) -cut_freedom = Engine.dimension(cut_ideal_ab_s) -println("Two points on a sphere, after cut: ", cut_freedom, " degrees of freedom") -if cut_freedom == 0 - coordring = base_ring(ideal_ab_s) - vbls = Variable.(symbols(coordring)) - cut_system = System([eqns_ab_s; cut], variables = vbls) - ##cut_result = HomotopyContinuation.solve(cut_system) - ##println("non-singular solutions:") - ##for soln in solutions(cut_result) - ## display(soln) - ##end - ##println("singular solutions:") - ##for sing in singular(cut_result) - ## display(sing.solution) - ##end - - # test a random witness set - max_slope = 2 - binom = Binomial(2max_slope, 1/2) - Random.seed!(6071) - samples = [] - for _ in 1:3 - cut_matrix = rand(binom, freedom, length(gens(coordring))) .- max_slope - ##cut_matrix = [ - ## 1 1 1 1 0 1 1 0 1 1 0; - ## 1 2 1 2 0 1 1 0 1 1 0; - ## 1 1 0 1 0 1 2 0 2 0 0 - ##] - sph_z_ind = indexin([sph.coords[5] for sph in ctx.spheres], gens(coordring)) - cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)] - println("sphere z variables: ", vbls[sph_z_ind]) - display(cut_matrix) - 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 +coordring = base_ring(ideal_ab_s) +vbls = Variable.(symbols(coordring)) +##cut_system = CompiledSystem(System([eqns_ab_s; cut], variables = vbls)) +##cut_result = HomotopyContinuation.solve(cut_system) +##println("non-singular solutions:") +##for soln in solutions(cut_result) +## display(soln) +##end +##println("singular solutions:") +##for sing in singular(cut_result) +## display(sing.solution) +##end + +# test a random witness set +system = CompiledSystem(System(eqns_ab_s, variables = vbls)) +max_slope = 2 +binom = Binomial(2max_slope, 1/2) +Random.seed!(6071) +samples = [] +for _ in 1:3 + cut_matrix = rand(binom, freedom, length(gens(coordring))) .- max_slope + ##cut_matrix = [ + ## 1 1 1 1 0 1 1 0 1 1 0; + ## 1 2 1 2 0 1 1 0 1 1 0; + ## 1 1 0 1 0 1 2 0 2 0 0 + ##] + sph_z_ind = indexin([sph.coords[5] for sph in ctx.spheres], gens(coordring)) + cut_offset = [sum(cf[sph_z_ind]) for cf in eachrow(cut_matrix)] + println("sphere z variables: ", vbls[sph_z_ind]) + display(cut_matrix) + display(cut_offset) + cut_subspace = LinearSubspace(cut_matrix, cut_offset) + wtns = witness_set(system, 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 = 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))") + 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 \ No newline at end of file