Diophantine Polynomial 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 x inside is called a bound occurrence; any other occurrence of x is called a free occurrence. A variable x is called a free variable of a polynomial system if the system contains a free occurrence of x. 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 notebook 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.

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 Developer`HermiteNormalForm. 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.

This system has four variables and two independent equations, hence the result is expressed in terms of two integer parameters.

In[1]:=

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), backsubstitutes 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 nonnegative integer parameters, and the number of new parameters may be larger than the number of variables.

This system has three variables; however, to express the solution set, we need eight nonnegative integer parameters.

In[2]:=

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.

Here Reduce finds integer solutions of a sparse polynomial of degree 100,000.

In[3]:=

The free term of this polynomial has 609,152 digits and cannot be easily factored.

In[6]:=

In[7]:=

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.

Here Reduce finds integer solutions of and selects the ones that satisfy the inequality .

In[8]:=

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 DiscreteSolutionBound. 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 (3) 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.

Hyperbolic type equations with square determinants

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

Here is a binary quadratic equation with = 9.

In[9]:=

Hyperbolic type equations with non-square determinants

If and is not an integer, then the equation (3) is a Pell-type equation. Methods for solving such equations have been developed since the 18th 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.

A Pell equation is an equation of the form , where D is not a square. Solutions to Pell equations with small coefficients can be quite complicated.

In[10]:=

Here is the solution of a Pell-type equation with = 5.

In[11]:=

Even though the solutions are expressed using nonrational numbers, they are in fact integers, as they should be.

In[12]:=

Reduce can solve systems consisting of a Pell-type equation and inequalities giving simple bounds on variables.

In[13]:=

Parabolic type equations

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

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

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

If , then the equation (4) implies

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

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

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

Here Reduce finds integer solutions of a quadratic equation of the parabolic type.

In[14]:=

Elliptic type equations

If , the solutions of equation (3) 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].

Here Reduce finds positive integer solutions of a quadratic equation of the elliptic type. There are more than possible positive integer values of , so the obvious algorithm would not be practical for this ellipse.

In[15]:=

Thue Equations

A Thue equation is a Diophantine equation of the form

where is an irreducible homogenous form of degree ≥ 3.

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.

This finds integer solutions of a cubic Thue equation.

In[16]:=

If we give Reduce a bound on the size of solutions, it can solve the equation much faster.

In[17]:=

Here Reduce finds the only solution of a degree 15 Thue equation with at most a 100-digit coordinate. Without the bound on the solution size, Reduce did not finish in 1000 seconds.

In[18]:=

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

This equation has no solutions modulo 2.

In[19]:=

Here FindInstance finds a small solution using the modular sieve.

In[20]:=

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, which 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.

Here the first equation has four integer solutions for x and y. For each of the solutions, the second equation becomes a univariate equation in z.

In[21]:=

Here the first equation is a Thue equation with one solution. After replacing x and y with the values computed from the first equation, the second equation becomes a Pell equation.

In[22]:=

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

substitutes the solutions into the equation and the inequalities present in the original system, and attempts to solve the so obtained systems for the integer parameters present in the solutions of the system (8).

This system reduces to solving a congruence and a Pell equation.

In[23]:=

This system reduces to solving a system of two congruences and a quadratic Diophantine equation of the parabolic type for each family of congruence solutions.

In[24]:=

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 (9). 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 (9).

This solves a sum of squares equation in three variables.

In[25]:=

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.

This cubic equation reduces to 12 pairs of quadratic and linear equations.

In[26]:=

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 (10) is solved using the method for solving systems with more than one equation. Mathematica recognizes three cases when the second part of the system (10) 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 (10) 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 (10) 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 (10) can be selected from a finite number of solution candidates.

Additionally, Mathematica uses the following heuristic to detect cases when the system (10) has no solutions. If there is an integer , such that is always divisible by m, and is never divisible by m, then the system (10) has no solutions. Candidates for m are found by computing the GCD of the values of f 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 SieveMaxPoints.

This reduces to (10) with .

In[27]:=

This reduces to (10) with .

In[28]:=

This reduces to the case of system (10).

In[29]:=

Here Reduce detects that the equation has no solutions, because is always divisible by 5, and is never divisible by 5.

In[30]:=

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, we enumerate all possible integer values of the first coordinate, then for each value of the first coordinate, we 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. By default this bound is equal to 10. It can be changed using the system option DiscreteSolutionBound. For systems for which the real solution set is unbounded or bounded but large, the solution is represented implicitly by returning 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].

Here the real solution set is bounded, but Reduce gives some cylinders in an implicit form. This is because some of the intervals bounding y contain more than 10 integers.

In[31]:=

Increasing the value of the system option DiscreteSolutionBound allows Reduce to find all integer solutions explicitly.

In[32]:=

This resets DiscreteSolutionBound back to the default value.

In[34]:=

Here the modular sieve method shows that there are no solutions in . After adding the corresponding inequalities, Reduce shows that this equation has no solutions.

In[35]:=

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 (11) 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.

This transforms to a system (11) with no real solutions.

In[36]:=

Here the system (11) obtained after three recursive transformations has a reducible nonconstant part.

In[37]:=

Options

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

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 f[1], f[2], ...

By default, the new parameters generated by Reduce are named C[1], C[2], ...

In[38]:=

The option GeneratedParameters allows users to customize the parameter names.

In[39]:=

ReduceOptions Group of System Options

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

DiscreteSolutionBound

The value of the system option DiscreteSolutionBound 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 n. The default value of the option is 10.

There are 10 integers in the real interval . Reduce writes them out explicitly.

In[40]:=

There are 11 integers in the real interval . Reduce represents them implicitly.

In[41]:=

This increases the DiscreteSolutionBound to 11.

In[42]:=

Now Reduce represents the solutions explicitly.

In[43]:=

This resets DiscreteSolutionBound back to the default value.

In[44]:=

SieveMaxPoints

The system option SieveMaxPoints 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 default value of the option is 10000.

With the default setting of SieveMaxPoints, FindInstance is unable to find a solution for this equation.

In[45]:=

Increasing the number of SieveMaxPoints to one million allows FindInstance to find a solution.

In[46]:=

This resets SieveMaxPoints to the default value.

In[48]:=

References

[1] C. Goldbach, Letter to L. Euler, June 7, 1742. http://www.mathstat.dal.ca/~joerg/pic/g-letter.jpg or http://www.informatik.uni-giessen.de/staff/richstein/pic/g-letter-zoomed.jpg.

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

[3] E. Contejean and H. Devie, "An efficient incremental algorithm for solving systems of linear Diophantine equations," Information and Computation, 113, 1994, pp. 143-172.

[4] F. Cucker, P. Koiran, and S. Smale, "A Polynomial Time Algorithm for Diophantine Equations in One Variable," J. Symb. Comput., 27 (1), 1999, pp. 21-30.

[5] A. Strzebonski, "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] L. E. Dickson, History of the Theory of Numbers, New York: Chelsea, 1952.

[7] T. Nagell, Introduction to Number Theory, New York: Wiley, 1951.

[8] K. Hardy, J. B. Muskat and K. S. Williams, "A deterministic algorithm for solving in coprime integers and ," Math. Comput., 55 (1990), pp. 327-343.

[9] N. Smart, The Algorithmic Resolution of Diophantine Equations, London Mathematical Society Student Text, 41, Cambridge University Press, 1998.

[10] S. Wagon, "The Magic of Imaginary Factoring", Mathematica in Education and Research, 5:1 (1996), pp. 43-47.