# Diophantine Polynomial Systems

Introduction | Multivariate Nonlinear Systems |

Linear Systems | Options |

Univariate Systems | References |

Bivariate Systems |

## Introduction

A *Diophantine polynomial system* is an expression constructed with polynomial equations and inequalities

combined using logical connectives and quantifiers

where the variables represent integer quantities.

An occurrence of a variable inside or is called a *bound occurrence*; any other occurrence of is called a *free occurrence.* A variable is called a *free variable* of a polynomial system if the system contains a free occurrence of . A Diophantine polynomial system is *quantifier-free* if it contains no quantifiers. A *decision problem* is a system with all variables existentially quantified, that is, a system of the form

where are all variables in . The decision problem (1) is equivalent to True or False, depending on whether the quantifier-free system of polynomial equations and inequalities has integer solutions.

An example of a Diophantine polynomial system is

Goldbach's conjecture [1], formulated in 1742 and still unproven, asserts that system (2) is equivalent to True. This suggests that *Mathematica* may not be able to solve arbitrary Diophantine polynomial systems. In fact, Matiyasevich's solution of Hilbert's tenth problem [2] shows that no algorithm can be constructed that would solve arbitrary Diophantine polynomial systems, not even quantifier-free systems or decision problems. Nevertheless, *Mathematica *functions Reduce, Resolve, and FindInstance are able to solve several reasonably large classes of Diophantine systems. This tutorial describes these classes of systems and methods used by *Mathematica* to solve them. The methods are presented in the order in which they are used. The tutorial also covers options affecting the methods used and how they operate.

## Linear Systems

### Systems of Linear Equations

Conjunctions of linear Diophantine equations are solvable for an arbitrary number of variables. *Mathematica* uses a method based on the computation Hermite normal form of matrices, available in *Mathematica* directly as HermiteDecomposition. The result may contain new unrestricted integer parameters. If the equations are independent, the number of parameters is equal to the difference between the number of variables and the number of equations.

In[1]:= |

Out[1]= |

### Frobenius Equations

A *Frobenius equation* is an equation of the form

where are positive integers, is an integer, and the coordinates of solutions are required to be non-negative integers.

For finding solution instances of Frobenius equations *Mathematica* uses a fast algorithm based on the computation of the critical tree in the Frobenius graph [11]. The algorithm applies when the smallest of does not exceed the value of the system option (the default is 1,000,000). Otherwise the more general methods for solving bounded linear systems are used. Functions FrobeniusSolve and FrobeniusNumber provide specialized functionality for solving Frobenius equations and computing Frobenius numbers.

In[2]:= |

Out[2]= |

### Bounded Systems of Linear Equations and Inequalities

If a real solution set of a system of linear equations and inequalities is a bounded polyhedron, the system has finitely many integer solutions. To find the solutions, *Mathematica* uses the following procedure.

You may assume the system has the form , where is a integer matrix, is a length integer vector, is an integer matrix, and is a length integer vector. First, the method for solving systems of linear equations is used to find an integer vector such that and a integer matrix *N *whose rows generate the null space of . The integer solution set of is equal to . Put and . The integer solution set of is equal to , where is the integer solution set of . To improve efficiency of finding the set , *Mathematica* simplifies using LatticeReduce, obtaining . Note that if the columns of are linearly dependent, has no solutions (otherwise it would have infinitely many solutions, which contradicts the assumptions). Hence you may assume that has linearly independent columns and so has columns. Put . Lattice reduction techniques are also used to find a small vector in the lattice . Let be such that . The set can be computed from the set of all such that using the formula .

To find the set a simple recursive algorithm can be used. The algorithm finds the bounds on the first variable using LinearProgramming and, for each integer value between the bounds, calls itself recursively with the first variable set to . This algorithm is used when the system option is set to False. With the default setting True a hybrid algorithm combining the recursive algorithm and a branch-and-bound-type algorithm is used. At each level of the recursion, the recursion is continued for the "middle" values of the first variable (defined as a projection of the set of points contained in the real solution set together with a unit cube) while the remaining parts of the real solution set are searched for integer solutions using the branch-and-bound type algorithm. FindInstance finds the single element of it needs using a branch-and-bound type algorithm.

There are two system options, and , that allow you to control the exact algorithm used. In some cases changing the values of these options may greatly improve the performance of Reduce.

In[3]:= |

Out[3]= |

*Mathematica* enumerates the solutions explicitly only if the number of integer solutions of the system does not exceed the maximum of the power of the value of the system option , where is the dimension of the solution lattice of the equations, and the second element of the value of the system option .

In[4]:= |

Out[4]= |

In[5]:= |

In[6]:= |

Out[6]= |

In[7]:= |

### Arbitrary Systems of Linear Equations and Inequalities

Quantifier-free systems of linear Diophantine equations and inequalities are solvable for an arbitrary number of variables. The system is written in the *disjunctive normal form*, that is, as a disjunction of conjunctions, and each conjunction is solved separately. If a conjunction contains only equations, the method for solving systems of linear equations is used. If the difference between the number of variables and the number of equations is at most one, *Mathematica* solves the equations using the method for solving systems of linear equations, and if the solution contains at most one free parameter (which is true in the generic case), back substitutes the solution into the inequalities to determine inequality restrictions for the parameter. For all other quantifier-free systems of linear Diophantine equations and inequalities *Mathematica* uses the algorithm described in [3], with some linear-programming-based improvements for handling bounded variables. The result may contain new non-negative integer parameters, and the number of new parameters may be larger than the number of variables.

In[8]:= |

Out[8]= |

## Univariate Systems

### Univariate Equations

To find integer solutions of univariate equations *Mathematica* uses a variant of the algorithm given in [4] with improvements described in [5]. The algorithm can find integer solutions of polynomials of much higher degrees than can be handled by real root isolation algorithms and with much larger free terms than can be handled by integer factorization algorithms.

In[9]:= |

Out[11]= |

In[12]:= |

Out[12]= |

In[13]:= |

Out[13]= |

### Systems of Univariate Equations and Inequalities

Systems of univariate Diophantine equations and inequalities are written in the disjunctive normal form, and each conjunction is solved separately. If a conjunction contains an equation, the method for solving univariate equations is used, and the solutions satisfying the remaining equations and inequalities are selected.

In[14]:= |

Out[14]= |

Conjunctions containing only inequalities are solved over the reals. Integer solutions in the resulting real intervals are given explicitly if their number in the given interval does not exceed the value of the system option . The default value of the option is 10. For intervals containing more integer solutions, the solutions are represented implicitly.

## Bivariate Systems

### Quadratic Equations

*Mathematica* can solve arbitrary quadratic Diophantine equations in two variables. The general form of such an equation is

If , where and are linear polynomials, the equation (1) is equivalent to , and methods for solving linear Diophantine equations are used. For irreducible polynomials , the algorithms used and the form of the result depend on the determinant of the quadratic form. The algorithms may use integer factorization and hence the correctness of the results depends on the correctness of the probabilistic primality test used by PrimeQ.

#### Hyperbolic-Type Equations with Square Determinants

If and is an integer, then , where and are linear polynomials and . In this case, the equation (1) is equivalent to the disjunction of linear systems , for all divisors of . Each of the linear systems has one solution over the rationals, hence the equation (1) has a finite number of integer solutions.

In[15]:= |

Out[15]= |

#### Hyperbolic-Type Equations with Nonsquare Determinants

If and is not an integer, then the equation (1) is a Pell-type equation. Methods for solving such equations have been developed since the 18 century and can be constructed based on [6] and [7] (though these books do not contain a complete description of the algorithm). The solution set is empty or infinite, parametrized by an integer parameter appearing in the exponent.

In[16]:= |

Out[16]= |

In[17]:= |

Out[17]= |

In[18]:= |

Out[18]= |

In[19]:= |

Out[19]= |

#### Parabolic-Type Equations

If , set , , and . Since , and are nonzero integers, and . Then

Set and . Then the equation (1) is equivalent to

Suppose . If the equation (1) had integer solutions, would have integer solutions in , and so would be a product of two linear polynomials. Since here is irreducible, the equation (1) has no integer solutions.

If , then the equation (2) implies

If the modular equation (3) has no solutions in , the equation (1) has no integer solutions. (If , the modular equation (3) has one solution, .) Otherwise , for some solution of the modular equation (3). Replacing in the equation (2) and solving the resulting linear equation for gives

Note that since satisfies the modular equation (3), the division in the last term of (4) gives an integer result. Since and , . Taking the equation (4) and the fact that into account gives

Therefore, the full solution of the equation (1) of parabolic type consists of one-parameter solution families given by equations (4) and (5) for each solution of the modular equation (3), for which is an integer.

In[20]:= |

Out[20]= |

#### Elliptic-Type Equations

If , the solutions of equation (1) are integer points on an ellipse. Since an ellipse is a bounded set, the number of solutions must be finite. An obvious algorithm for finding integer points would be to compute the solutions for for each of the finite number of possible integer values of . This, however, would be prohibitively slow for larger ellipses. *Mathematica* uses a faster algorithm described in [8].

In[21]:= |

Out[21]= |

### Thue Equations

A* Thue *equation is a Diophantine equation of the form

where is an irreducible homogeneous form of degree

The number of solutions of Thue equations is always finite. *Mathematica* can in principle solve arbitrary Thue equations, though the time necessary to find the solutions lengthens very fast with degree and coefficient size. The hardest part of the algorithm is computing a bound on the size of solutions. *Mathematica* uses an algorithm based on the Baker-Wustholz theorem to find the bound [9]. If the input contains inequalities that provide a reasonable size bound on solutions, *Mathematica* can compute the solutions much faster.

In[22]:= |

Out[22]= |

In[23]:= |

Out[23]= |

In[24]:= |

Out[24]= |

## Multivariate Nonlinear Systems

### Systems Solvable with the Modular Sieve Method

*Mathematica* uses a variant of the modular sieve method (see e.g. [9]). The method may prove that a system has no solutions in integers modulo an integer , and therefore, it has no integer solutions. Otherwise, it may find a solution with small integer coordinates or prove that the system has no integer solutions with all coordinates between and . The number of candidate solution points that the sieve method is allowed to test is controlled by the system option .

In[25]:= |

Out[25]= |

In[26]:= |

Out[26]= |

### Systems with More than One Equation

If a Diophantine polynomial system contains more than one equation, *Mathematica* uses GroebnerBasis in an attempt to reduce the problem to a sequence of simpler problems.

#### Systems Solvable by Recursion over Finitely Many Partial Solutions

*Mathematica* attempts to solve an element of the Gröbner basis that depends on the minimal number of the initial variables. If the number of solutions is finite, then for each solution *Mathematica* substitutes the computed values and attempts to solve the obtained system for the remaining variables.

In[27]:= |

Out[27]= |

In[28]:= |

Out[28]= |

#### Systems with Linear-Triangular Gröbner Bases

This heuristic applies to systems with Gröbner bases of the form

In this case, *Mathematica* solves the system of congruences

The solutions are represented using new integer parameters. These are substituted into the equation and the inequalities present in the original system, and *Mathematica* attempts to solve the so-obtained systems for the new parameters.

In[29]:= |

Out[29]= |

In[30]:= |

Out[30]= |

### Sums of Squares

*Mathematica* can solve equations of the form

using the algorithm described in [10]. For multivariate quadratic equations, *Mathematica* attempts to find an affine transformation that puts the equation in the form (2). A heuristic method based on CholeskyDecomposition is used for this purpose. However, the method may fail for some equations that can be represented in the form (2).

In[31]:= |

Out[31]= |

To find a single solution of (2) FindInstance uses an algorithm based on [12].

In[32]:= |

Out[32]= |

In[33]:= |

Out[33]= |

### Pythagorean Equation

*Mathematica* knows the solution to the Pythagorean equation

In[34]:= |

Out[34]= |

For quadratic equations in three variables, *Mathematica* attempts to find a transformation of the form

transforming the equation to the Pythagorean equation.

In[35]:= |

Out[35]= |

### Equations with Reducible Nonconstant Parts

If the sum of nonconstant terms in an equation factors, *Mathematica* uses the formula

to reduce the equation to a disjunction of pairs of equations with lower degrees. Note that this reduction depends on the ability to find all divisors of , hence the correctness of the results depends on the correctness of the probabilistic primality test used by PrimeQ.

In[36]:= |

Out[36]= |

### Equations with a Linear Variable

*Mathematica* attempts to solve Diophantine systems of the form

where is a conjunction of inequalities or True, by reducing them to

The first part of the system (3) is solved using the method for solving systems with more than one equation. *Mathematica* recognizes three cases when the second part of the system (3) is solvable. If , the solution is given by and by the restrictions on obtained by solving the inequalities . Nonlinear systems of inequalities are solved using CylindricalDecomposition. If for an integer constant , the solution of the second part of the system (3) is given by and by the restrictions on obtained by solving the congruence and then solving the inequalities for each solution of the congruence. If is nonconstant, *Mathematica* can solve the second part of the system (3) if . Since *Mathematica* factors all equations at the preprocessing stage, and can be assumed to be relatively prime. Then

for an integer and polynomials and with integer coefficients and . If is an integer, then is an integer, and so or . Since , the last condition is satisfied only by a finite number of integers . Hence the solutions of the second part of the system (3) can be selected from a finite number of solution candidates.

Additionally, *Mathematica* uses the following heuristic to detect cases when the system (3) has no solutions. If there is an integer , such that is always divisible by , and is never divisible by , then the system (3) has no solutions. Candidates for are found by computing the GCD of the values of at several points.

The last two methods use exhaustive search over finite sets of points. The allowed number of search points is controlled by the system option .

In[37]:= |

Out[37]= |

In[38]:= |

Out[38]= |

In[39]:= |

Out[39]= |

In[40]:= |

Out[40]= |

### Systems with Empty or Bounded Real Solution Sets

If a Diophantine polynomial system is not solved by any other methods, *Mathematica* solves the system over the reals using the Cylindrical Algebraic Decomposition (CAD) algorithm. If the system has no real solutions, then clearly it has no integer solutions. If the real solution set is bounded, then the number of integer solutions is finite. In principle, all the integer solutions can be found in this case from a cylindrical decomposition. Namely, for each cylinder, you enumerate all possible integer values of the first coordinate, then for each value of the first coordinate, you enumerate all possible integer values of the second coordinate, and so on. However, for large bounded solution sets this method could lead to a huge number of points to try. Therefore, *Mathematica* has a bound on the number of explicitly enumerated integer solutions in a real interval. For systems for which the real solution set is unbounded or bounded but large, the solution is represented implicitly by returning the CAD and a condition that all variables are integers. Note that for multivariate systems such an implicit representation may not even be enough to tell whether integer solutions exist. This should be expected, given Matiyasevich's solution of Hilbert's tenth problem [2].

In[41]:= |

Out[41]= |

In[42]:= |

Out[42]= |

In[43]:= |

Out[44]= |

In[45]:= |

In[46]:= |

Out[46]= |

### Equations of the Form

*Mathematica* attempts to solve Diophantine systems of the form

where is a conjunction of inequalities or True, by transforming them to

The resulting system (4) may, or may not, be easier to solve. Systems exist for which this transformation could be applied recursively arbitrarily many times; therefore, *Mathematica* uses a recursion bound to ensure the heuristic terminates.

In[47]:= |

Out[47]= |

In[48]:= |

Out[48]= |

### Systems Solvable by Exhaustive Search

For systems containing explicit lower and upper bounds on all variables, *Mathematica* uses exhaustive search to find solutions. The bounds of the search are specified by the value of the system option . The option value should be a pair of integers (the default is ). If the number of integer points within the bounds does not exceed the first integer, the exhaustive search is used instead of any solution methods other than univariate polynomial solving. Otherwise, if the number of integer points within the bounds does not exceed the second integer, the exhaustive search is performed after all other methods fail.

In[49]:= |

Out[49]= |

## Options

The *Mathematica* functions for solving Diophantine polynomial systems have a number of options that control the way they operate. This tutorial gives a summary of these options.

option name | default value | |

GeneratedParameters | C | specifies how the new parameters generated to represent solutions should be named |

Reduce options affecting the behavior for Diophantine polynomial systems.

### GeneratedParameters

To represent infinite solutions of some Diophantine systems, Reduce needs to introduce new integer parameters. The names of the new parameters are specified by the option GeneratedParameters. With GeneratedParameters->f, the new parameters are named , , ... .

In[50]:= |

Out[50]= |

In[51]:= |

Out[51]= |

### ReduceOptions Group of System Options

Here are the system options from the group that may affect the behavior of Reduce, Resolve, and FindInstance for Diophantine polynomial systems. The options can be set with

SetSystemOptions["ReduceOptions"->{"option name"->value}].

option name | default value | |

"BranchLinearDiophantine" | True | whether Reduce should use a branch-and-bound-type algorithm to compute solutions of bounded systems of linear Diophantine inequalities |

"DiscreteSolutionBound" | 10 | the bound on the number of explicitly enumerated integer solutions in a real interval |

"ExhaustiveSearchMaxPoints" | {1000,10000} | the maximal number of integer points within variable bounds for which the exhaustive search is used before and after all other solution methods |

"ImplicitIntegerSolutions" | Automatic | whether Reduce should be allowed to return real solutions with Element conditions when it cannot find explicit integer solutions |

"LatticeReduceDiophantine" | True | whether LatticeReduce should be used to preprocess bounded systems of linear Diophantine inequalities |

"MaxFrobeniusGraph" | 1000000 | the maximal size of the smallest coefficient in a Frobenius equation for which FindInstance computes the critical tree in the Frobenius graph |

"SieveMaxPoints" | {10000,1000000} | the maximal number of points at which the modular sieve method evaluates the system before and after all other solution methods |

group options affecting the behavior of Reduce, Resolve, and FindInstance for Diophantine polynomial systems.

#### BranchLinearDiophantine

The value of the system option specifies which variant of the algorithm should be used in the final stage of solving bounded linear systems. With the default setting "BranchLinearDiophantine"->False, a simple recursive enumeration is used. With "BranchLinearDiophantine"->True, Reduce uses a hybrid method combining a branch-and-bound-type algorithm and a simple recursive enumeration.

In[52]:= |

Out[54]= |

In[55]:= |

Out[56]= |

In[57]:= |

Out[66]= |

In[67]:= |

Out[68]= |

#### DiscreteSolutionBound

The value of the system option specifies whether integer solutions in a real interval should be enumerated explicitly or represented implicitly as . With DiscreteSolutionBound->n, the integer solutions in the given real interval are enumerated explicitly if their number does not exceed . The default value of the option is 10.

In[69]:= |

Out[69]= |

In[70]:= |

Out[70]= |

In[71]:= |

In[72]:= |

Out[72]= |

The value of also affects the solving of bounded linear systems and other Diophantine systems with finitely many solutions. When you need an explicit enumeration of a large number of solutions in a bounded region, set to Infinity.

In[73]:= |

Out[75]= |

In[76]:= |

Out[77]= |

In[78]:= |

#### ExhaustiveSearchMaxPoints

The system option specifies the maximal number of search points used by the exhaustive search method. The option value should be a pair of integers (the default is ). If the number of integer points within the bounds does not exceed the first integer, the exhaustive search is used instead of any solution methods other than univariate polynomial solving. Otherwise, if the number of integer points within the bounds does not exceed the second integer, the exhaustive search is performed after all other methods fail.

In[79]:= |

Out[79]= |

In[80]:= |

In[81]:= |

Out[81]= |

In[82]:= |

Out[82]= |

In[83]:= |

Out[84]= |

In[85]:= |

Out[86]= |

In[87]:= |

Out[87]= |

In[88]:= |

#### ImplicitIntegerSolutions

The value of the system option specifies whether integer solutions in a real interval should be enumerated explicitly or represented implicitly as . With DiscreteSolutionBound->n, the integer solutions in the given real interval are enumerated explicitly if their number does not exceed . The default value of the option is 10.

In[89]:= |

Out[89]= |

In[90]:= |

Out[91]= |

In[92]:= |

Out[93]= |

In[94]:= |

Out[95]= |

In[96]:= |

#### LatticeReduceDiophantine

The value of the system option specifies whether LatticeReduce should be used to preprocess systems of bounded linear inequalities. The use of LatticeReduce is important for systems of inequalities describing polyhedra whose projections on some nonaxial lines are much smaller than their projections on the axes. However, there are systems for which LatticeReduce, instead of simplifying the problem, makes it significantly harder.

In[97]:= |

Out[98]= |

In[99]:= |

In[100]:= |

Out[100]= |

*,*which bound solutions to a reasonably small size polyhedron, combined with a set of relatively complicated inequalities . For such systems, using LatticeReduce tends to increase the timing.

In[101]:= |

Out[107]= |

In[108]:= |

Out[109]= |

In[110]:= |

#### MaxFrobeniusGraph

The system option specifies the maximal size of the smallest coefficient in a Frobenius equation for which FindInstance uses an algorithm based on the computation of the critical tree in the Frobenius graph [11]. Otherwise, the more general methods for solving bounded linear systems are used. Unlike the general method for solving bounded linear systems, the method based on the computation of the Frobenius graph depends very little on the number of variables, hence it is the faster choice for equations with many variables. On the other hand, the method requires storing a graph of the size of the smallest coefficient, so for large coefficients it may run out of memory.

In[1]:= |

Out[4]= |

In[5]:= |

Out[5]= |

In[6]:= |

In[7]:= |

Out[7]= |

In[8]:= |

Out[8]= |

In[9]:= |

#### SieveMaxPoints

The system option specifies the maximal number of search points used by the modular sieve method and by searches used in solving equations with a linear variable. The value of the option is a pair . Search involving up to points is tried before some of the slower solving methods. For decision problems, another search involving up to points is tried after all other methods fail. The default value of the option is .

In[10]:= |

Out[10]= |

In[11]:= |

Out[11]= |

In[12]:= |

Out[13]= |

In[14]:= |

Out[14]= |

In[15]:= |

Out[15]= |

In[16]:= |

Out[17]= |

In[18]:= |

## References

[1] Goldbach, C. Letter to L. Euler, June 7, 1742. http://mathworld.wolfram.com/GoldbachConjecture.html. [online version]

[2] Matiyasevich, Yu. V. "The Diophantineness of Enumerable Sets (Russian)" *Dokl. Akad. Nauk SSSR* 191 (1970): 279-282. English translation: *Soviet Math. Dokl.* 11 (1970): 354-358.

[3] Contejean, E. and H. Devie. "An Efficient Incremental Algorithm for Solving Systems of Linear Diophantine Equations." *Information and Computation* 113 (1994): 143-172.

[4] Cucker, F., P. Koiran, and S. Smale. "A Polynomial Time Algorithm for Diophantine Equations in One Variable." *Journal of Symbolic Computation* 27, no. 1 (1999): 21-30.

[5] Strzebonski, A. "An Improved Algorithm for Diophantine Equations in One Variable." Paper presented at the ACA 2002 Session on Symbolic-Numerical Methods in Computational Science, Volos, Greece. Notebook with the conference talk available at members.wolfram.com/adams.

[6] Dickson, L. E. *History of the Theory of Numbers*. Chelsea, 1952.

[7] Nagell, T. *Introduction to Number Theory.* Wiley, 1951.

[8] Hardy, K., J. B. Muskat and K. S. Williams. "A Deterministic Algorithm for Solving in Coprime Integers and ." *Mathematics of Computation* 55 (1990): 327-343.

[9] Smart, N. *The Algorithmic Resolution of Diophantine Equations*. Cambridge University Press, 1998.

[10] Bressoud, D. M. and S. Wagon. *A Course in Computational Number Theory*. Key College Publishing, 2000.

[11] Beihoffer, D. E., J. Hendry, A. Nijenhuis, and S. Wagon. "Faster Algorithms for Frobenius Numbers." to appear in *The Electronic Journal of Combinatorics.*

[12] Rabin, M. O. and J. O. Shallit. "Randomized Algorithms in Number Theory."* Communications on Pure and Applied Mathematics* 39 (1986): 239-256.