# Showing that an integer solution does not exist

I am given a quadratic form (non-homogeneous) in 12 variables and I want to show that it does not have an integer solution.

The 'FindInstance' and 'Solve' command both give {} while 'Reduce' returns 'false'. So that is good news. However, I want to use the fact in a paper, and simply relying on a mathematica result is not the good way to go in algebraic geometry.

My next step would have been to consider the equation modulo different prime numbers and have mathematica check for solutions there. As soon as I find a prime number such that the equation does not have a solution modulo that prime, I could start to think about an argument there. However, starting from p=5, the program refuses to run. Can an equation be so complicated that the program does not even try to solve it, or what could be the issue here?

I'm using the following command:

 FindInstance[
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) == -2, {a, b, c, d, e, f, g, h, i, j, k, l},
Modulus -> 5]


Additonally: apart from Simplify and FullSimplify, does mathematica have special commands for quadratic forms to transform them into a standard form?

• Hi Louis. Do you already know that it does not have integer solutions, and you are just trying to prove it? Or are you trying to determine if it does not have integer solutions? – MikeY Apr 24 '17 at 13:11
• Hi MikeY. I do not know for sure if it does not have integer solutions. However, the fact that 'Reduce' delivers 'False' as an answer should prove exactly that modulo algorithmic errors in mathematica. Also, I have heuristical reasons to believe that there are no integer solutions. – Louis Apr 24 '17 at 13:13
• @Louis Is it possible that this system has only complex solutions? Try NSolve. – zhk Apr 24 '17 at 13:51
• @MMM it has very simple rational solutions (for example c=-1/2, e=-1/2, everything else=0 is a solution) – Louis Apr 24 '17 at 13:53
• @Louis I thought that idea of finding one prime modulo where there is no solution implies no integer solutions was only for linear diophantine equations. Yours are non linear. – bobbym Apr 24 '17 at 15:18

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

• Thank you very much for your extremely enlightening and helpful answer. I actually had not considered the Cholesky Decomposition so far - a computational tool I last saw 5 years ago in a class of my bachelor studies. I think the following argument might work in the specific situation: we take the Cholesky Decomposition (as proposed by you). Now we observe that assuming there exists a solution for our equation, the $k$-component has absolut value at least one (simply by considering the problem modulo 4). But the $k$-component of cd.ivec only depends on the $k$ and $l$ components of ivec, and – Louis Apr 26 '17 at 9:44
• from the entries in cd it should be clear that for absolute value of the $k$-component of ivec at least 1, we either have that the $k$ or $l$ component of cd.ivec already has absolute value bigger than $\sqrt{2}$ which should give a short proof. – Louis Apr 26 '17 at 9:46
• Of course, the rest of your answer is still extremely helpful to me, because I will have to compute more similar problems and at least I now have an idea about how mathematicas (or sages) algorithms could work. – Louis Apr 26 '17 at 9:48
• I realized after posting that I had both an error and a missed optimization. The error is that not all variables should be constrained to integer values; this is only true for the ones that denote multiples of the unit vectors. The missed optimization is along the lines of one of these comments, to the effect that indivudual components of cd.ivec may be given have size bounds. I edited with these considerations taken into account. – Daniel Lichtblau Apr 26 '17 at 15:23
• The observation abs(k)>=1 is certainly useful. But it is not clear to me how it will force the 11th or 12th entry of cd.ivec to exceed sqrt(2) in magnitude. – Daniel Lichtblau Apr 26 '17 at 15:35