Try out the Gram matrix approach

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

@@ -0,0 +1,30 @@
+using LinearAlgebra
+using AbstractAlgebra
+
+F, (a, b, c) = rational_function_field(Generic.Rationals{BigInt}(), ["a", "b", "c"])
+M = matrix_space(F, 5, 5)
+
+# 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;
+]))
+
+r, p, L, U = lu(gram)
+solution = transpose(p * L)
+mform = U * inv(transpose(L))
+
+concrete = [evaluate(entry, [0, 1, -3//4]) for entry in solution]
+
+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_concrete = std_basis * concrete

engine-proto/gram-test/gram-test.sage

@@ -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