engine-proto/gram-test/gram-test.jl

@@ -1,23 +1,56 @@
 using LinearAlgebra
 using AbstractAlgebra
 
-F, (a, b, c) = rational_function_field(Generic.Rationals{BigInt}(), ["a", "b", "c"])
-M = matrix_space(F, 5, 5)
+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 = M(F.([
-  -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;
+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(2-n_pos); -ones(3-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(gram)
-solution = transpose(p * L)
-mform = U * inv(transpose(L))
+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)
 
-concrete = [evaluate(entry, [0, 1, -3//4]) for entry in 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;
@@ -27,4 +60,18 @@ std_basis = [
   0 0 1 0  0
 ]
 std_solution = M(F.(std_basis)) * solution
-std_concrete = std_basis * concrete
+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)