Here is something of a brute-force approach.
poly = c*(-4*c + 2*d) + d*(2*c - 4*d + 2*e) +
e*(2*d - 4*e + 2*f) + (-2*b + 2*f + 2*g - 4*h)*h +
f*(2*e - 4*f + 2*h) + g*(-4*g + 2*h) +
j*(126*b - 28*i - 196*j + 182*k) +
i*(56*a + 154*b - 140*i - 28*j - 63*k - 560*l) +
k*(42*a - 28*b - 63*i + 182*j - 266*k - 420*l) + (238*a + 560*b -
560*i - 420*k - 2380*l)*l +
a*(-24*a - 56*b + 56*i + 42*k + 238*l) +
b*(-56*a - 220*b - 2*h + 154*i + 126*j - 28*k + 560*l);
vars = Variables[poly];
Create the corresponding symmetric quadratic form.
qform = Normal[CoefficientArrays[poly, vars][[3]]];
qform2 = (qform + Transpose[qform])/2;
Check this:
Expand[vars.qform2.vars - poly]
(* Out[161]= 0 *)
It is negative definite (easy to check: all eigenvalues are negative). Find the Cholesky decomposition of its negative.
cd = CholeskyDecomposition[-qform2];
Check this:
Max[Abs[Transpose[cd].cd + qform2]]
(* Out[163]= 0 *)
Integer solutions of the form ivec will satisfy ivec.qform2.ivec==-2. This in turn means we have a fundamental equation:
(cd.ivec) . (cd.ivec) == 2
If it can be shown that no integer vector ivec can satisfy this, then that does the job. We'll show how to (fairly tightly) bound the search space, after which one can do an exhaustive search.
Note that any integer vector that works is of course an integer combination of the 12 unit vectors (our quadratic form is 12x12 because we have 12 variables in the polynomial). For a vector of the form cd.ivec, if it has 1-norm larger than or equal to 5 then a standard argument shows the 2-norm (Euclidean length) must be at least 5/12*sqrt(12) which is in turn slightly larger than sqrt(2) (so the dot product of this vector with itself would exceed 2). (For fixed 1-norm, the 2-norm is minimized where all components are equal). We can use this to bound the search space. Actually we can tighten to 2 sqrt(6) plus a small increment to allow for approximate arithmetic.
Another perhaps stronger constraint is that no individual component can exceed sqrt(2) in magnitude.
The optimization step uses mixed linear integer programming, with integer variables for the multiples of the unit vectors. We use a constraint that a certain one (in a loop) be at least 1. As noted above we also constraint that the 1-norm not exceed 2 sqrt(6)+epsilon and that individual components not exceed sqrt(2) in size. We use the fact that negating ivec gives a vector cd.ivec of all opposite signs (thus not changing the fundamental equation) so we only need compute upper bounds. I do the computation at machine precision so a strictly symbolic proof would require redoing the ILP in exact arithmetic (more than I'll do but by no means impossible).
We use cvars below to form the integer vector, absvars to denote absolute values of components of cd.cvars, impose constraints for sizes, and iteratively force each component of cvars to be greater-equal to one. Note that absvars need not be integer valued.
cvars = Array[cc, Length[vars]];
absvars = Array[aa, Length[vars]];
cvec = cd.cvars;
c1 = Thread[cvec <= absvars];
c2 = Thread[-cvec <= absvars];
c3 = Map[0 <= # <= Sqrt[2] &, absvars];
c4 = Total[absvars] <= 2*Sqrt[6] + 1/1000;
c5 = Element[cvars, Integers];
bounds = Table[c6 = cc[j] >= 1;
NMaximize[{cc[j], Join[c1, c2, c3, {c4, c5, c6}]},
Join[cvars, absvars]], {j, Length[cvars]}];
bounds[[All, 1]]
(* Out[220]= {1., 2., 3., 4., 2., 2., 4., 1., 2., 1., 1., 1.} *)
This give a large, but computationally tractable, search space (<100M cases, by my count). One can perhaps improve matters by breaking into more special cases though. Here is our vector whose norm-square needs to be 2.
cvec
(* Out[244]= {2 cc[1] - cc[2], Sqrt[3] cc[2] - (2 cc[3])/Sqrt[3],
2 Sqrt[2/3] cc[3] - Sqrt[3/2] cc[4],
Sqrt[5/2] cc[4] - 2 Sqrt[2/5] cc[7],
2 Sqrt[55] cc[5] + cc[7]/Sqrt[55] - (63 cc[8])/Sqrt[55] -
7 Sqrt[11/5] cc[9] + (14 cc[10])/Sqrt[55] + (28 cc[11])/Sqrt[55] -
56 Sqrt[5/11] cc[12], 2 cc[6] - cc[7],
2 Sqrt[19/55] cc[7] + (63 cc[8])/(2 Sqrt[1045]) +
7/2 Sqrt[11/95] cc[9] - (7 cc[10])/Sqrt[1045] - (14 cc[11])/Sqrt[
1045] + 28 Sqrt[5/209] cc[12],
7/2 Sqrt[953/95] cc[8] - (3331 cc[9])/(2 Sqrt[90535]) - (
4499 cc[10])/Sqrt[90535] + (882 cc[11])/Sqrt[90535] -
1764 Sqrt[5/18107] cc[12],
2 Sqrt[34/953] cc[9] + (89 cc[10])/(2 Sqrt[32402]) -
14 Sqrt[2/16201] cc[11] + 140 Sqrt[2/16201] cc[12],
1/2 Sqrt[5271/34] cc[10] - 26 Sqrt[14/12801] cc[11] +
260 Sqrt[14/12801] cc[12],
2 Sqrt[38/753] cc[11] - (7 cc[12])/Sqrt[28614], Sqrt[7/38] cc[12]} *)
The first and sixth components are fairly convenient in form.
cvec[[{1, 6}]]
(* Out[246]= {2 cc[1] - cc[2], 2 cc[6] - cc[7]} *)
Neither can exceed 1 in size so one can create 3x3=9 separate cases from:
{2*cc[1] - 1 <= cc[2] <= 2*cc[1] + 1,
2*cc[6] - 1 <= cc[7] <= 2*cc[6] + 1}
Not sure whether this will help substantially. The original formulation I had would benefit but this improved one might not.
--- edit ---
We can get further improvement by simply noting that any sum of two elements in cd.ivec cannot exceed 2 (because the 2-norm would be minimized when each was 1, and that 2-norm squared would then be the target value of 2. This will tighten the search space.
cvars = Array[cc, Length[vars]];
absvars = Array[aa, Length[vars]];
cvec = cd.cvars;
c1 = Thread[cvec <= absvars];
c2 = Thread[-cvec <= absvars];
c3 = Map[0 <= # <= Sqrt[2] &, absvars];
c4 = Total[absvars] <= 2*Sqrt[6] + 1/1000;
c5 = Element[cvars, Integers];
c7 = Map[0 <= # <= 2 &,
Union[Flatten[Outer[Plus, absvars, absvars]]] /. 2*_ :> Nothing];
bounds = Table[c6 = cc[j] >= 1;
NMaximize[{cc[j], Join[c1, c2, c3, {c4, c5, c6}, c7]},
Join[cvars, absvars]], {j, Length[cvars]}];
bounds[[All, 1]]
(* Out[231]= {1., 1., 2., 3., 2., 2., 3., 1., 2., 1., 1., 1.} *)
This brings down the size of the search space.
Round[Times @@ Map[2*# + 1 &, bounds[[All, 1]]]]
(* Out[232]= 22325625 *)
One can do similarly with triples of the absolute value variables, restricting sums to be no larger than 3 sqrt(2/3). This gives a healthy improvement, bringing the bounds to {1., 1., 1., 1., 1., 1., 1., 2., 1., 1., 1., 1.}, for a search space of less than 10^6 integer combinations for the original variables.
--- end edit ---
NSolve. $\endgroup$