Some statistics about Diophantine Equations

Andrej Bauer

November 3, 2006

We use Mathematica to compute statistics about solvability of small Diophantine equations in natural numbers.

Generation of terms

We are going to generate terms in the language of Peano arithmetic with zero, successor, addition and multiplication.  First we define the successor function:

ClearAll[s]
s[n_] := n + 1

Functions terms generates terms of a given complexity with given variables (this will always include all terms of smaller complexity, up to equality). Our notion of complexity is the depth of the expression tree. At each step of generation we put the terms into canonical form and throw out duplicates.

ClearAll[combine, terms,  
      termsUpTo]

combine[f_, lst_] := 
Flatten[
        Function[t1, 
            Function[t2, f[t1, t2]] /@ Select[lst, OrderedQ[{t1, #1}] &]] /@ 
          lst, 1]

terms[0, vars_] := terms[0, vars] = 
        Prepend[vars, 0]
terms[n_,  
        vars_] := terms[n, 
          vars] = 
Union [
Map[
            Expand@s[#] &,  terms[n - 1, vars]], 
combine[
            Expand@Plus[#1, #2] &,  
            terms[n - 1, vars]], 
combine[
            Expand@Times[#1, #2] &,  
            terms[n - 1, vars]]
]

Generation of equations

Given a list of terms lst, we generate all equations (Subscript[e, 1]) = (Subscript[e, 2]), where we avoid taking both (Subscript[e, 1]) = (Subscript[e, 2]) and (Subscript[e, 2]) = (Subscript[e, 1]).

ClearAll[equations]
equations[n_, 
        vars_] := 
With[{lst = terms[n, vars]},  
        Union @ 
          Flatten[Function[e1,  
                Function[e2,  e1 == e2] /@  
                  Select[lst, OrderedQ[{e1, #}] &]]  /@  lst, 
            1]]

We can use Mathematica function FullSimplify to rewrite equations in equivalent form. Many equations are reduced to the same form this way. In the following example the original equations are in the left column and the simplified ones in the right column.

{#, 
            FullSimplify[#, 
              Assumptions -> {x ∈ Integers, 
                  y ∈ Integers, x ≥ 0, 
                  y ≥ 0}]} & /@ 
        Take[equations[2, {x, y}], {1, 4000, 200}] // TableForm

Function reducedEquations generates all reduced equations of given depth in given variables.

ClearAll[reducedEquations]
reducedEquations[n_, vars_] := 
      Union@FullSimplify[equations[n, vars], 
          Assumptions -> {And @@ (# ∈ 
                          Integers && # ≥ 0 & /@ 
                    vars)}]

The number of equations of depth 2 in two variables is almost halved by this procedure:

equations[2, {x, y}] // Length

reducedEquations[2, {x, y}] // Length

Diophantine Equations

Mathematica can solve a Diophantine equation in natural numbers with the FindInstance function:

FindInstance[{(((1 + x))^2) + (((1 + y))^2) == (((1 + z))^2) && 
          x ≥ 0 && y ≥ 0 && 
          z ≥ 0}, {x, y, z}, Integers]

It returns an empty list when there is no solution:

FindInstance[{(((2  x + 1))^2) == 4  y + 3 && 
          x ≥ 0 && y ≥ 0}, {x, y}, 
      Integers]

It cannot deal with arbitrary equations, of course:

FindInstance[{(((1 + x))^7) + (((1 + y))^7) == (((1 + z))^7) && 
          x ≥ 0 && y ≥ 0 && 
          z ≥ 0}, {x, y, z}, Integers]

We define a function which returns True, False or Maybe depending on whether a given equation has a solution:

ClearAll[hasSolution]
Off[FindInstance::nsmet]
hasSolution[eq_, vars_] := 
      Switch[
FindInstance[
          eq && And @@ (# ≥ 0 & /@ vars), vars, 
          Integers], 
{},  
        False, 
_List,  
        True, 
_FindInstance,  
        Print["Cannot decide: ",  eq]; Maybe]

Examples of usage:

hasSolution[(((1 + x))^2) + (((1 + y))^2) == (((1 + z))^2), {x, y, z}]

hasSolution[(((2  x + 1))^2) == 4  y + 3, {x, y}]

hasSolution[(((1 + x))^7) + (((1 + y))^7) == (((1 + z))^7), {x, y, z}]

Now we can write a function which collects statistics about a set of equations:

ClearAll[stats]
stats[eqs_, vars_] := 
      With[{r = hasSolution[#, vars] & /@ 
              eqs}, 
{{True, Count[r, True]}, {False, 
            Count[r, False]}, {Maybe, Count[r, Maybe]}}]

As an example, let us look at statistic for equations with term depth 1 in one variable:

terms[1, {x}]

equations[1, {x}] // ColumnForm

Of the above 16 equations, 11 have a solution and 5 do not:

stats[equations[1, {x}], {x}] // TableForm

The same example with reduced equations:

reducedEquations[1, {x}] // ColumnForm

Now we have only 7 equations left, two of which do not have a solution:

stats[reducedEquations[1, {x}], {x}] // TableForm

Statistics

Depth 1

One variable:

stats[equations[1, {x}], {x}] // TableForm

Two variables:

stats[equations[1, {x, y}], {x, y}] // TableForm

Three variables:

stats[equations[1, {x, y, z}], {x, y, z}] // TableForm

Four variables:

stats[equations[1, {x, y, z, w}], {x, y, z, w}] // TableForm

Five variables:

stats[equations[1, {x, y, z, w, v}], {x, y, z, w, v}] // 
      TableForm

Depth 2

One variable:

stats[equations[2, {x}], {x}] // TableForm

Two variables:

stats[equations[2, {x, y}], {x, y}] // TableForm

Three variables:

stats[equations[2, {x, y, z}], {x, y, z}] // TableForm

Depth 3

Depth 3 is quite large. We can only handle the one variable case:

stats[equations[3, {x}], {x}] // TableForm