Compare commits
38 Commits
Author | SHA1 | Date |
---|---|---|
Aaron Fenyes | ef33b8ee10 | |
Aaron Fenyes | 717e5a6200 | |
Aaron Fenyes | 16826cf07c | |
Aaron Fenyes | 3170a933e4 | |
Aaron Fenyes | f2000e5731 | |
Aaron Fenyes | ba365174d3 | |
Aaron Fenyes | ae5db0f9ea | |
Aaron Fenyes | 8d8bc9162c | |
Aaron Fenyes | 291d5c8ff6 | |
Aaron Fenyes | e41bcc7e13 | |
Aaron Fenyes | 31d5e7e864 | |
Aaron Fenyes | a450f701fb | |
Aaron Fenyes | 6cf07dc6a1 | |
Aaron Fenyes | 1f173708eb | |
Aaron Fenyes | 6f18d4efcc | |
Aaron Fenyes | 621c4c5776 | |
Aaron Fenyes | b3b7c2026d | |
Aaron Fenyes | af1d31f6e6 | |
Aaron Fenyes | 8e33987f59 | |
Aaron Fenyes | 06872a04af | |
Aaron Fenyes | becefe0c47 | |
Aaron Fenyes | 34358a8728 | |
Aaron Fenyes | 95c0ff14b2 | |
Aaron Fenyes | f97090c997 | |
Aaron Fenyes | 45aaaafc8f | |
Aaron Fenyes | 43cbf8a3a0 | |
Aaron Fenyes | 21f09c4a4d | |
Aaron Fenyes | a3f3f6a31b | |
Aaron Fenyes | 65d23fb667 | |
Aaron Fenyes | 4e02ee16fc | |
Aaron Fenyes | 6349f298ae | |
Aaron Fenyes | 0731c7aac1 | |
Aaron Fenyes | 59a527af43 | |
Aaron Fenyes | c29000d912 | |
Aaron Fenyes | 86dbd9ea45 | |
Aaron Fenyes | 463a3b21e1 | |
Aaron Fenyes | 4d5aa3b327 | |
Aaron Fenyes | b864cf7866 |
|
@ -0,0 +1,203 @@
|
|||
module Algebraic
|
||||
|
||||
export
|
||||
codimension, dimension,
|
||||
Construction, realize,
|
||||
Element, Point, Sphere,
|
||||
Relation, LiesOn, AlignsWithBy, mprod
|
||||
|
||||
import Subscripts
|
||||
using LinearAlgebra
|
||||
using AbstractAlgebra
|
||||
using Groebner
|
||||
using ...HittingSet
|
||||
|
||||
# --- commutative algebra ---
|
||||
|
||||
# as of version 0.36.6, AbstractAlgebra only supports ideals in multivariate
|
||||
# polynomial rings when the coefficients are integers. we use Groebner to extend
|
||||
# support to rationals and to finite fields of prime order
|
||||
Generic.reduce_gens(I::Generic.Ideal{U}) where {T <: FieldElement, U <: MPolyRingElem{T}} =
|
||||
Generic.Ideal{U}(base_ring(I), groebner(gens(I)))
|
||||
|
||||
function codimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}}
|
||||
leading = [exponent_vector(f, 1) for f in gens(I)]
|
||||
targets = [Set(findall(.!iszero.(exp_vec))) for exp_vec in leading]
|
||||
length(HittingSet.solve(HittingSetProblem(targets), maxdepth))
|
||||
end
|
||||
|
||||
dimension(I::Generic.Ideal{U}, maxdepth = Inf) where {T <: RingElement, U <: MPolyRingElem{T}} =
|
||||
length(gens(base_ring(I))) - codimension(I, maxdepth)
|
||||
|
||||
# --- primitve elements ---
|
||||
|
||||
abstract type Element{T} end
|
||||
|
||||
mutable struct Point{T} <: Element{T}
|
||||
coords::Vector{MPolyRingElem{T}}
|
||||
vec::Union{Vector{MPolyRingElem{T}}, Nothing}
|
||||
rel::Nothing
|
||||
|
||||
## [to do] constructor argument never needed?
|
||||
Point{T}(
|
||||
coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[],
|
||||
vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing
|
||||
) where T = new(coords, vec, nothing)
|
||||
end
|
||||
|
||||
function buildvec!(pt::Point)
|
||||
coordring = parent(pt.coords[1])
|
||||
pt.vec = [one(coordring), dot(pt.coords, pt.coords), pt.coords...]
|
||||
end
|
||||
|
||||
mutable struct Sphere{T} <: Element{T}
|
||||
coords::Vector{MPolyRingElem{T}}
|
||||
vec::Union{Vector{MPolyRingElem{T}}, Nothing}
|
||||
rel::Union{MPolyRingElem{T}, Nothing}
|
||||
|
||||
## [to do] constructor argument never needed?
|
||||
Sphere{T}(
|
||||
coords::Vector{MPolyRingElem{T}} = MPolyRingElem{T}[],
|
||||
vec::Union{Vector{MPolyRingElem{T}}, Nothing} = nothing,
|
||||
rel::Union{MPolyRingElem{T}, Nothing} = nothing
|
||||
) where T = new(coords, vec, rel)
|
||||
end
|
||||
|
||||
function buildvec!(sph::Sphere)
|
||||
coordring = parent(sph.coords[1])
|
||||
sph.vec = sph.coords
|
||||
sph.rel = mprod(sph.coords, sph.coords) + one(coordring)
|
||||
end
|
||||
|
||||
const coordnames = IdDict{Symbol, Vector{Union{Symbol, Nothing}}}(
|
||||
nameof(Point) => [nothing, nothing, :xₚ, :yₚ, :zₚ],
|
||||
nameof(Sphere) => [:rₛ, :sₛ, :xₛ, :yₛ, :zₛ]
|
||||
)
|
||||
|
||||
coordname(elt::Element, index) = coordnames[nameof(typeof(elt))][index]
|
||||
|
||||
function pushcoordname!(coordnamelist, indexed_elt::Tuple{Any, Element}, coordindex)
|
||||
eltindex, elt = indexed_elt
|
||||
name = coordname(elt, coordindex)
|
||||
if !isnothing(name)
|
||||
subscript = Subscripts.sub(string(eltindex))
|
||||
push!(coordnamelist, Symbol(name, subscript))
|
||||
end
|
||||
end
|
||||
|
||||
function takecoord!(coordlist, indexed_elt::Tuple{Any, Element}, coordindex)
|
||||
elt = indexed_elt[2]
|
||||
if !isnothing(coordname(elt, coordindex))
|
||||
push!(elt.coords, popfirst!(coordlist))
|
||||
end
|
||||
end
|
||||
|
||||
# --- primitive relations ---
|
||||
|
||||
abstract type Relation{T} end
|
||||
|
||||
mprod(v, w) = (v[1]*w[2] + w[1]*v[2]) / 2 - dot(v[3:end], w[3:end])
|
||||
|
||||
# elements: point, sphere
|
||||
struct LiesOn{T} <: Relation{T}
|
||||
elements::Vector{Element{T}}
|
||||
|
||||
LiesOn{T}(pt::Point{T}, sph::Sphere{T}) where T = new{T}([pt, sph])
|
||||
end
|
||||
|
||||
equation(rel::LiesOn) = mprod(rel.elements[1].vec, rel.elements[2].vec)
|
||||
|
||||
# elements: sphere, sphere
|
||||
struct AlignsWithBy{T} <: Relation{T}
|
||||
elements::Vector{Element{T}}
|
||||
cos_angle::T
|
||||
|
||||
AlignsWithBy{T}(sph1::Sphere{T}, sph2::Sphere{T}, cos_angle::T) where T = new{T}([sph1, sph2], cos_angle)
|
||||
end
|
||||
|
||||
equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) - rel.cos_angle
|
||||
|
||||
# --- constructions ---
|
||||
|
||||
mutable struct Construction{T}
|
||||
points::Vector{Point{T}}
|
||||
spheres::Vector{Sphere{T}}
|
||||
relations::Vector{Relation{T}}
|
||||
|
||||
function Construction{T}(; elements = Vector{Element{T}}(), relations = Vector{Relation{T}}()) where T
|
||||
allelements = union(elements, (rel.elements for rel in relations)...)
|
||||
new{T}(
|
||||
filter(elt -> isa(elt, Point), allelements),
|
||||
filter(elt -> isa(elt, Sphere), allelements),
|
||||
relations
|
||||
)
|
||||
end
|
||||
end
|
||||
|
||||
function Base.push!(ctx::Construction{T}, elt::Point{T}) where T
|
||||
push!(ctx.points, elt)
|
||||
end
|
||||
|
||||
function Base.push!(ctx::Construction{T}, elt::Sphere{T}) where T
|
||||
push!(ctx.spheres, elt)
|
||||
end
|
||||
|
||||
function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T
|
||||
push!(ctx.relations, rel)
|
||||
for elt in rel.elements
|
||||
push!(ctx, elt)
|
||||
end
|
||||
end
|
||||
|
||||
function realize(ctx::Construction{T}) where T
|
||||
# collect coordinate names
|
||||
coordnamelist = Symbol[]
|
||||
eltenum = enumerate(Iterators.flatten((ctx.spheres, ctx.points)))
|
||||
for coordindex in 1:5
|
||||
for indexed_elt in eltenum
|
||||
pushcoordname!(coordnamelist, indexed_elt, coordindex)
|
||||
end
|
||||
end
|
||||
|
||||
# construct coordinate ring
|
||||
coordring, coordqueue = polynomial_ring(parent_type(T)(), coordnamelist, ordering = :degrevlex)
|
||||
|
||||
# retrieve coordinates
|
||||
for (_, elt) in eltenum
|
||||
empty!(elt.coords)
|
||||
end
|
||||
for coordindex in 1:5
|
||||
for indexed_elt in eltenum
|
||||
takecoord!(coordqueue, indexed_elt, coordindex)
|
||||
end
|
||||
end
|
||||
|
||||
# construct coordinate vectors
|
||||
for (_, elt) in eltenum
|
||||
buildvec!(elt)
|
||||
end
|
||||
|
||||
# turn relations into equations
|
||||
eqns = vcat(
|
||||
equation.(ctx.relations),
|
||||
[elt.rel for (_, elt) in eltenum if !isnothing(elt.rel)]
|
||||
)
|
||||
|
||||
# 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
|
||||
|
||||
(Generic.Ideal(coordring, eqns), eqns)
|
||||
end
|
||||
|
||||
end
|
|
@ -0,0 +1,53 @@
|
|||
module Numerical
|
||||
|
||||
using Random: default_rng
|
||||
using LinearAlgebra
|
||||
using AbstractAlgebra
|
||||
using HomotopyContinuation:
|
||||
Variable, Expression, AbstractSystem, System, LinearSubspace,
|
||||
nvariables, isreal, witness_set, results
|
||||
import GLMakie
|
||||
using ..Algebraic
|
||||
|
||||
# --- polynomial conversion ---
|
||||
|
||||
# hat tip Sascha Timme
|
||||
# https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl/issues/520#issuecomment-1317681521
|
||||
function Base.convert(::Type{Expression}, f::MPolyRingElem)
|
||||
variables = Variable.(symbols(parent(f)))
|
||||
f_data = zip(coefficients(f), exponent_vectors(f))
|
||||
sum(cf * prod(variables .^ exp_vec) for (cf, exp_vec) in f_data)
|
||||
end
|
||||
|
||||
# create a ModelKit.System from an ideal in a multivariate polynomial ring. the
|
||||
# variable ordering is taken from the polynomial ring
|
||||
function System(I::Generic.Ideal)
|
||||
eqns = Expression.(gens(I))
|
||||
variables = Variable.(symbols(base_ring(I)))
|
||||
System(eqns, variables = variables)
|
||||
end
|
||||
|
||||
# --- sampling ---
|
||||
|
||||
function real_samples(F::AbstractSystem, dim; rng = default_rng())
|
||||
# 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)...
|
||||
))
|
||||
cut = LinearSubspace(normals, fill(0., dim))
|
||||
filter(isreal, results(witness_set(F, cut, seed = 0x1974abba)))
|
||||
end
|
||||
|
||||
AbstractAlgebra.evaluate(pt::Point, vals::Vector{<:RingElement}) =
|
||||
GLMakie.Point3f([evaluate(u, vals) for u in pt.coords])
|
||||
|
||||
function AbstractAlgebra.evaluate(sph::Sphere, vals::Vector{<:RingElement})
|
||||
radius = 1 / evaluate(sph.coords[1], vals)
|
||||
center = radius * [evaluate(u, vals) for u in sph.coords[3:end]]
|
||||
GLMakie.Sphere(GLMakie.Point3f(center), radius)
|
||||
end
|
||||
|
||||
end
|
|
@ -0,0 +1,76 @@
|
|||
include("HittingSet.jl")
|
||||
|
||||
module Engine
|
||||
|
||||
include("Engine.Algebraic.jl")
|
||||
include("Engine.Numerical.jl")
|
||||
|
||||
using .Algebraic
|
||||
using .Numerical
|
||||
|
||||
export Construction, mprod, codimension, dimension
|
||||
|
||||
end
|
||||
|
||||
# ~~~ sandbox setup ~~~
|
||||
|
||||
using Random
|
||||
using Distributions
|
||||
using LinearAlgebra
|
||||
using AbstractAlgebra
|
||||
using HomotopyContinuation
|
||||
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)
|
||||
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")
|
||||
|
||||
# --- test rational cut ---
|
||||
|
||||
coordring = base_ring(ideal_tan_sph)
|
||||
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")
|
||||
|
||||
# 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
|
|
@ -0,0 +1,111 @@
|
|||
module HittingSet
|
||||
|
||||
export HittingSetProblem, solve
|
||||
|
||||
HittingSetProblem{T} = Pair{Set{T}, Vector{Pair{T, Set{Set{T}}}}}
|
||||
|
||||
# `targets` should be a collection of Set objects
|
||||
function HittingSetProblem(targets, chosen = Set())
|
||||
wholeset = union(targets...)
|
||||
T = eltype(wholeset)
|
||||
unsorted_moves = [
|
||||
elt => Set(filter(s -> elt ∉ s, targets))
|
||||
for elt in wholeset
|
||||
]
|
||||
moves = sort(unsorted_moves, by = pair -> length(pair.second))
|
||||
Set{T}(chosen) => moves
|
||||
end
|
||||
|
||||
function Base.display(problem::HittingSetProblem{T}) where T
|
||||
println("HittingSetProblem{$T}")
|
||||
|
||||
chosen = problem.first
|
||||
println(" {", join(string.(chosen), ", "), "}")
|
||||
|
||||
moves = problem.second
|
||||
for (choice, missed) in moves
|
||||
println(" | ", choice)
|
||||
for s in missed
|
||||
println(" | | {", join(string.(s), ", "), "}")
|
||||
end
|
||||
end
|
||||
println()
|
||||
end
|
||||
|
||||
function solve(pblm::HittingSetProblem{T}, maxdepth = Inf) where T
|
||||
problems = Dict(pblm)
|
||||
while length(first(problems).first) < maxdepth
|
||||
subproblems = typeof(problems)()
|
||||
for (chosen, moves) in problems
|
||||
if isempty(moves)
|
||||
return chosen
|
||||
else
|
||||
for (choice, missed) in moves
|
||||
to_be_chosen = union(chosen, Set([choice]))
|
||||
if isempty(missed)
|
||||
return to_be_chosen
|
||||
elseif !haskey(subproblems, to_be_chosen)
|
||||
push!(subproblems, HittingSetProblem(missed, to_be_chosen))
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
problems = subproblems
|
||||
end
|
||||
problems
|
||||
end
|
||||
|
||||
function test(n = 1)
|
||||
T = [Int64, Int64, Symbol, Symbol][n]
|
||||
targets = Set{T}.([
|
||||
[
|
||||
[1, 3, 5],
|
||||
[2, 3, 4],
|
||||
[1, 4],
|
||||
[2, 3, 4, 5],
|
||||
[4, 5]
|
||||
],
|
||||
# example from Amit Chakrabarti's graduate-level algorithms class (CS 105)
|
||||
# notes by Valika K. Wan and Khanh Do Ba, Winter 2005
|
||||
# https://www.cs.dartmouth.edu/~ac/Teach/CS105-Winter05/
|
||||
[
|
||||
[1, 3], [1, 4], [1, 5],
|
||||
[1, 3], [1, 2, 4], [1, 2, 5],
|
||||
[4, 3], [ 2, 4], [ 2, 5],
|
||||
[6, 3], [6, 4], [ 5]
|
||||
],
|
||||
[
|
||||
[:w, :x, :y],
|
||||
[:x, :y, :z],
|
||||
[:w, :z],
|
||||
[:x, :y]
|
||||
],
|
||||
# Wikipedia showcases this as an example of a problem where the greedy
|
||||
# algorithm performs especially poorly
|
||||
[
|
||||
[:a, :x, :t1],
|
||||
[:a, :y, :t2],
|
||||
[:a, :y, :t3],
|
||||
[:a, :z, :t4],
|
||||
[:a, :z, :t5],
|
||||
[:a, :z, :t6],
|
||||
[:a, :z, :t7],
|
||||
[:b, :x, :t8],
|
||||
[:b, :y, :t9],
|
||||
[:b, :y, :t10],
|
||||
[:b, :z, :t11],
|
||||
[:b, :z, :t12],
|
||||
[:b, :z, :t13],
|
||||
[:b, :z, :t14]
|
||||
]
|
||||
][n])
|
||||
problem = HittingSetProblem(targets)
|
||||
if isa(problem, HittingSetProblem{T})
|
||||
println("Correct type")
|
||||
else
|
||||
println("Wrong type: ", typeof(problem))
|
||||
end
|
||||
problem
|
||||
end
|
||||
|
||||
end
|
|
@ -0,0 +1,85 @@
|
|||
using LinearAlgebra
|
||||
using AbstractAlgebra
|
||||
|
||||
function printgood(msg)
|
||||
printstyled("✓", color = :green)
|
||||
println(" ", msg)
|
||||
end
|
||||
|
||||
function printbad(msg)
|
||||
printstyled("✗", color = :red)
|
||||
println(" ", msg)
|
||||
end
|
||||
|
||||
F, gens = rational_function_field(Generic.Rationals{BigInt}(), ["a₁", "a₂", "b₁", "b₂", "c₁", "c₂"])
|
||||
a = gens[1:2]
|
||||
b = gens[3:4]
|
||||
c = gens[5:6]
|
||||
|
||||
# three mutually tangent spheres which are all perpendicular to the x, y plane
|
||||
gram = [
|
||||
-1 1 1;
|
||||
1 -1 1;
|
||||
1 1 -1
|
||||
]
|
||||
|
||||
eig = eigen(gram)
|
||||
n_pos = count(eig.values .> 0.5)
|
||||
n_neg = count(eig.values .< -0.5)
|
||||
if n_pos + n_neg == size(gram, 1)
|
||||
printgood("Non-degenerate subspace")
|
||||
else
|
||||
printbad("Degenerate subspace")
|
||||
end
|
||||
sig_rem = Int64[ones(1-n_pos); -ones(4-n_neg)]
|
||||
unk = hcat(a, b, c)
|
||||
M = matrix_space(F, 5, 5)
|
||||
big_gram = M(F.([
|
||||
diagm(sig_rem) unk;
|
||||
transpose(unk) gram
|
||||
]))
|
||||
|
||||
r, p, L, U = lu(big_gram)
|
||||
if isone(p)
|
||||
printgood("Found a solution")
|
||||
else
|
||||
printbad("Didn't find a solution")
|
||||
end
|
||||
solution = transpose(L)
|
||||
mform = U * inv(solution)
|
||||
|
||||
vals = [0, 0, 0, 1, 0, -3//4]
|
||||
solution_ex = [evaluate(entry, vals) for entry in solution]
|
||||
mform_ex = [evaluate(entry, vals) for entry in mform]
|
||||
|
||||
std_basis = [
|
||||
0 0 0 1 1;
|
||||
0 0 0 1 -1;
|
||||
1 0 0 0 0;
|
||||
0 1 0 0 0;
|
||||
0 0 1 0 0
|
||||
]
|
||||
std_solution = M(F.(std_basis)) * solution
|
||||
std_solution_ex = std_basis * solution_ex
|
||||
|
||||
println("Minkowski form:")
|
||||
display(mform_ex)
|
||||
|
||||
big_gram_recovered = transpose(solution_ex) * mform_ex * solution_ex
|
||||
valid = all(iszero.(
|
||||
[evaluate(entry, vals) for entry in big_gram] - big_gram_recovered
|
||||
))
|
||||
if valid
|
||||
printgood("Recovered Gram matrix:")
|
||||
else
|
||||
printbad("Didn't recover Gram matrix. Instead, got:")
|
||||
end
|
||||
display(big_gram_recovered)
|
||||
|
||||
# this should be a solution
|
||||
hand_solution = [0 0 1 0 0; 0 0 -1 2 2; 0 0 0 1 -1; 1 0 0 0 0; 0 1 0 0 0]
|
||||
unmix = Rational{Int64}[[1//2 1//2; 1//2 -1//2] zeros(Int64, 2, 3); zeros(Int64, 3, 2) Matrix{Int64}(I, 3, 3)]
|
||||
hand_solution_diag = unmix * hand_solution
|
||||
big_gram_hand_recovered = transpose(hand_solution_diag) * diagm([1; -ones(Int64, 4)]) * hand_solution_diag
|
||||
println("Gram matrix from hand-written solution:")
|
||||
display(big_gram_hand_recovered)
|
|
@ -0,0 +1,27 @@
|
|||
F = QQ['a', 'b', 'c'].fraction_field()
|
||||
a, b, c = F.gens()
|
||||
|
||||
# three mutually tangent spheres which are all perpendicular to the x, y plane
|
||||
gram = matrix([
|
||||
[-1, 0, 0, 0, 0],
|
||||
[0, -1, a, b, c],
|
||||
[0, a, -1, 1, 1],
|
||||
[0, b, 1, -1, 1],
|
||||
[0, c, 1, 1, -1]
|
||||
])
|
||||
|
||||
P, L, U = gram.LU()
|
||||
solution = (P * L).transpose()
|
||||
mform = U * L.transpose().inverse()
|
||||
|
||||
concrete = solution.subs({a: 0, b: 1, c: -3/4})
|
||||
|
||||
std_basis = matrix([
|
||||
[0, 0, 0, 1, 1],
|
||||
[0, 0, 0, 1, -1],
|
||||
[1, 0, 0, 0, 0],
|
||||
[0, 1, 0, 0, 0],
|
||||
[0, 0, 1, 0, 0]
|
||||
])
|
||||
std_solution = std_basis * solution
|
||||
std_concrete = std_basis * concrete
|
Loading…
Reference in New Issue