Mathematica quadric normalizer

quadrics.nb

K := 5 x^2 + 4 x y + 2 y^2 + 2 x - 4 y + 3 (* Input quadric *)

X := Variables[K]; n := Length[X]

(* Mathematica's Coefficient[3x+2xy,x] gets us 3+2y. We just want the free term, 3. *)
FreeTerm[p_] := p /. Map[# -> 0 &, X] (* p with all variables substituted by 0 *)
Co[m_] := FreeTerm[Coefficient[K, m]] (* So Co[x^2] is 5, Co[y] is -4... *)

A := Table[Co[X[[i]] * X[[j]]] / If[i==j, 1, 2], {i, 1, n}, {j, 1, n}]
B := Table[Co[X[[i]]], {i, 1, n}]
c := FreeTerm[K]

S := Transpose[Map[Normalize, Eigenvectors[A]]]
(* Kn is K, normalized to only have monomials x^2,y^2,... and x,y,... *)
Kn := Simplify[X.Transpose[S].A.S.X + B.S.X + c]

(* Now we eliminate the linear terms. *)
(* We turn ax^2+bx into a(x+b/2a)^2 - b^2/4a for each variable x. *)
Kf := Sum[Module[
   {x = X[[i]], a = Coefficient[Kn, X[[i]]^2], b = Coefficient[Kn, X[[i]]]},
   a (x+b/(2 a))^2 - b^2/(4a)],{i,1,n}] + c

Kf