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
_{x} or _{x} 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
quantifierfree if it contains no quantifiers. A
decision problem is a system with all variables existentially quantified, that is, a system of the form
where
x_{1}, ..., x_{n} are all variables in
. The decision problem (
1) is equivalent to
True or
False, depending on whether the quantifierfree system of polynomial equations and inequalities
(x_{1}, ..., x_{n}) 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 quantifierfree 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.
This system has four variables and two independent equations, hence the result is expressed in terms of two integer parameters.
Out[1]=  

Frobenius Equations
A
Frobenius equation is an equation of the form
where
a_{1}, ..., a_{n} are positive integers,
m is an integer, and the coordinates
x_{1}, ..., x_{n} of solutions are required to be nonnegative 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
a_{1}, ..., a_{n} does not exceed the value of the
MaxFrobeniusGraph 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.
This finds a solution of a Frobenius equation.
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
M_{eq}.x=b_{eq} M_{ineq}.x≥b_{ineq}, where
M_{eq} is a
kn integer matrix,
b_{eq} is a length
k integer vector,
M_{ineq} is an
ln integer matrix, and
b_{ineq} is a length
l integer vector. First,
the method for solving systems of linear equations is used to find an integer vector
s such that
M_{eq}.s=b_{eq} and a
pn integer matrix
N whose rows generate the nullspace of
M_{eq}.x=0. The integer solution set of
M_{eq}.x=b_{eq} is equal to
{s+i.N:i^{p}}. Put
M_{mult}=M_{ineq}.N^{T} and
b_{mult}=b_{ineq}M_{ineq}.s. The integer solution set of
M_{eq}.x=b_{eq} M_{ineq}.x≥b_{ineq} is equal to
{s+i.N:i}, where
is the integer solution set of
M_{mult}.i≥b_{mult}. To improve efficiency of finding the set
,
Mathematica simplifies
M_{mult}^{T} using
LatticeReduce, obtaining
M_{red}^{T}. Note that if the columns of
M_{mult} are linearly dependent,
M_{mult}.i≥b_{mult} has no solutions (otherwise it would have infinitely many solutions, which contradicts the assumptions). Hence you may assume that
M_{mult} has linearly independent columns and so
M_{red} has
p columns. Put
R= (M_{mult}^{T}.M_{mult})^{1} (M_{mult}^{T}M_{red}). Lattice reduction techniques are also used to find a small vector
b_{red} in the lattice
b_{mult}+M_{red}.v. Let
v_{0} be such that
b_{red}=b_{mult}+M_{red}.v_{0}. The set
can be computed from the set
_{red} of all
i^{p} such that
M_{red}.i≥b_{red} using the formula
={R. (iv_{0}):i_{red}}.
To find the set
_{red} a simple recursive algorithm can be used. The algorithm finds the bounds on the first variable using
LinearProgramming and, for each integer value
a_{1} between the bounds, calls itself recursively with the first variable set to
a_{1}. This algorithm is used when the system option
BranchLinearDiophantine is set to
False. With the default setting
True a hybrid algorithm combining the recursive algorithm and a branchandbound 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 branchandbound type algorithm.
FindInstance finds the single element of
_{red} it needs using a branchandbound type algorithm.
There are two system options,
BranchLinearDiophantine and
LatticeReduceDiophantine, 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.
This finds all integer points in a triangle.
Out[3]=  

Mathematica enumerates the solutions explicitly only if the number of integer solutions of the system does not exceed the maximum of the
p^{th} power of the value of the system option
DiscreteSolutionBound, where
p is the dimension of the solution lattice of the equations, and the second element of the value of the system option
ExhaustiveSearchMaxPoints.
Here Reduce does not give explicit solutions because their number would exceed the default limit of 10000.
Out[4]=  

This increases the value of the system option DiscreteSolutionBound to 1000. 
Since there are two variables and no equations, the limit on the number of solutions is now 1000^{2}, and Reduce can enumerate the solutions explicitly.
Out[6]=  

This resets DiscreteSolutionBound to the default value. 
Arbitrary Systems of Linear Equations and Inequalities
Quantifierfree 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 quantifierfree systems of linear Diophantine equations and inequalities
Mathematica uses the algorithm described in [
3], with some linearprogrammingbased 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, you need eight nonnegative integer parameters.
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.
Here Reduce finds integer solutions of a sparse polynomial of degree 100,000.
Out[11]=  

The free term of this polynomial has 609,152 digits and cannot be easily factored.
Out[12]=  
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.
Here Reduce finds integer solutions of x^{4}25 x^{2}144 and selects the ones that satisfy the inequality x^{100001}27x+5≥0.
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
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
(x, y)=_{1} (x, y)_{2} (x, y), where
_{1} (x, y) and
_{2} (x, y) are linear polynomials, the equation (
1) is equivalent to
_{1} (x, y)=0 _{2} (x, y)=0, and methods for solving linear Diophantine equations are used. For irreducible polynomials
(x, y), the algorithms used and the form of the result depend on the determinant
=b^{2}4 a c 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
>0 and
is an integer, then
(x, y)g=_{1} (x, y)_{2} (x, y), where
_{1} (x, y) and
_{2} (x, y) are linear polynomials and
g=c d^{2}+a e^{2}+b^{2}fb d e4 a c f. In this case, the equation (
1) is equivalent to the disjunction of linear systems
_{1} (x, y)= _{2} (x, y)=g/, for all divisors
of
g. Each of the linear systems has one solution over the rationals, hence the equation (
1) has a finite number of integer solutions.
Here is a binary quadratic equation with =9.
Out[15]=  

Hyperbolic Type Equations with Nonsquare Determinants
If
>0 and
is not an integer, then the equation (
1) is a Pelltype equation. Methods for solving such equations have been developed since the 18
^{th} 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 x^{2}D y^{2}=1, where D is not a square. Solutions to Pell equations with small coefficients can be quite complicated.
Out[16]=  

Here is the solution of a Pelltype equation with =5.
Out[17]=  

Even though the solutions are expressed using nonrational numbers, they are in fact integers, as they should be.
Out[18]=  

Reduce can solve systems consisting of a Pelltype equation and inequalities giving simple bounds on variables.
Out[19]=  

Parabolic Type Equations
If
=0, set
g=sign (a)gcd (a, c),
, and
. Since
b^{2}=4g^{2} (a/g) (c/g),
a_{1} and
c_{1} are nonzero integers, and
b=2g a_{1}c_{1}. Then
Set
m=c_{1}da_{1}e and
t=a_{1} x+c_{1} y. Then the equation (
1) is equivalent to
Suppose
m=0. If the equation (
1) had integer solutions,
a_{1}g t^{2}+d t+a_{1}f=0 would have integer solutions in
t, and so
(x, y) would be a product of two linear polynomials. Since here
(x, y) is irreducible, the equation (
1) has no integer solutions.
If
m≠0, then the equation (
2) implies
If the modular equation (
3) has no solutions in
t, the equation (
1) has no integer solutions. (If
m=1, the modular equation (
3) has one solution,
t=0.) Otherwise
t=u+k m, for some solution
0≤u<m of the modular equation (
3). Replacing
t→u+k m in the equation (
2) and solving the resulting linear equation for
y gives
Note that since
u satisfies the modular equation (
3), the division in the last term of (
4) gives an integer result. Since
t=a_{1} x+c_{1} y and
t=u+k m,
x= (u+k mc_{1}y)/a_{1}. Taking the equation (
4) and the fact that
m=c_{1}da_{1}e into account gives
Therefore, the full solution of the equation (
1) of parabolic type consists of oneparameter solution families given by equations (
4) and (
5) for each solution
u of the modular equation (
3), for which
(c_{1}g u^{2}+e u+c_{1}f)/m is an integer.
Here Reduce finds integer solutions of a quadratic equation of the parabolic type.
Out[20]=  

Elliptic Type Equations
If
<0, 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
y for each of the finite number of possible integer values of
x. 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 8×10^{18} possible positive integer values of x, so the obvious algorithm would not be practical for this ellipse.
Out[21]=  

Thue Equations
A
Thue equation is a Diophantine equation of the form
where
(x, y) 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 BakerWustholz 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.
Out[22]=  

If we give Reduce a bound on the size of solutions, it can solve the equation much faster.
Out[23]=  

Here Reduce finds the only solution of a degree 15 Thue equation with at most a 100digit x coordinate. Without the bound on the solution size, Reduce did not finish in 1000 seconds.
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
m, 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
b and
b. 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.
Out[25]=  

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.
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. The four univariate equations have a total of two integer solutions.
Out[27]=  

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.
Out[28]=  

Systems with LinearTriangular 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
g (Y)=0 and the inequalities present in the original system, and
Mathematica attempts to solve the soobtained systems for the new parameters.
This system reduces to solving a congruence and a Pell equation.
Out[29]=  

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.
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).
This solves a sumofsquares equation in three variables.
Out[31]=  

To find a single solution of (
2)
FindInstance uses an algorithm based on [
12].
This finds a decomposition of a 10000digit integer into a sum of seven squares. N is applied to make the printed result smaller.
Out[32]=  

This proves that the decomposition found is correct.
Out[33]=  

Pythagorean Equation
Mathematica knows the solution to the Pythagorean equation
This gives the general solution of the Pythagorean equation.
Out[34]=  

For quadratic equations in three variables,
Mathematica attempts to find a transformation of the form
transforming the equation to the Pythagorean equation.
This equation can be transformed to the Pythagorean equation.
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
c, hence the correctness of the results depends on the correctness of the probabilistic primality test used by
PrimeQ.
This cubic equation reduces to 12 pairs of quadratic and linear equations.
Out[36]=  

Equations with a Linear Variable
Mathematica attempts to solve Diophantine systems of the form
where
(x_{1}, ..., x_{n}, y) 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
f (x_{1}, ..., x_{n})1, the solution is given by
y=g (x_{1}, ..., x_{n}) and by the restrictions on
x_{1}, ..., x_{n} obtained by solving the inequalities
(x_{1}, ..., x_{n}, g (x_{1}, ..., x_{n})). Nonlinear systems of inequalities are solved using
CylindricalDecomposition. If
f (x_{1}, ..., x_{n})m for an integer constant
m≥2, the solution of the second part of the system (
3) is given by
y=g (x_{1}, ..., x_{n})/m and by the restrictions on
x_{1}, ..., x_{n} obtained by solving the congruence
g (x_{1}, ..., x_{n})0 mod m and then solving the inequalities
(x_{1}, ..., x_{n}, g (x_{1}, ..., x_{n})/m) for each solution of the congruence. If
f (x_{1}, ..., x_{n}) is nonconstant,
Mathematica can solve the second part of the system (
3) if
n=1. Since
Mathematica factors all equations at the preprocessing stage,
f (x_{1}) and
g (x_{1}) can be assumed to be relatively prime. Then
for an integer
d and polynomials
q (x_{1}) and
r (x_{1}) with integer coefficients and
deg (r)<deg (f). If
g (x_{1})/f (x_{1}) is an integer, then
r (x_{1})/f (x_{1}) is an integer, and so
r (x_{1})=0 or
r (x_{1})≥ f (x_{1}). Since
deg (r)<deg (f), the last condition is satisfied only by a finite number of integers
x_{1}. 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
m≥2, such that
f (x_{1}, ..., x_{n}) is always divisible by
m, and
g (x_{1}, ..., x_{n}) is never divisible by
m, then the system (
3) 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 ( 3) with f (x_{1}, ..., x_{n})1.
Out[37]=  

This reduces to ( 3) with f (x_{1}, ..., x_{n})3.
Out[38]=  

This reduces to the n=1 case of system ( 3).
Out[39]=  

Here Reduce detects that the equation has no solutions, because 9 x^{6} y^{3} z^{4}9 x^{2} y^{3} z^{8}5 y^{8} z^{9}10 is always divisible by 5, and 75 x^{4} y z^{4}+7 x^{8} y^{2} z^{4}9 z^{8}4 x^{6} y z^{8} is never divisible by 5.
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. 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 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].
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.
Out[41]=  

Increasing the value of the system option DiscreteSolutionBound allows Reduce to find all integer solutions explicitly.
Out[43]=  

This resets DiscreteSolutionBound to the default value. 
Here the modular sieve method shows that there are no solutions in (15, 15]^{3}. After adding inequalities to eliminate this cube, Reduce then recognizes that this equation has no solutions anywhere.
Out[45]=  

Equations of the Form x g (x, y, z_{1}, ..., z_{n})+y=c
Mathematica attempts to solve Diophantine systems of the form
where
(x, y, z_{1}, ..., z_{n}) 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.
This transforms to a system ( 4) with no real solutions.
Out[46]=  

Here the system ( 4) obtained after three recursive transformations has a reducible nonconstant part.
Out[47]=  

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
ExhaustiveSearchMaxPoints. The option value should be a pair of integers (the default is
{1000, 10000}). 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.
This transcendental Diophantine equation with bounded variable values is solved by exhaustive search.
Out[48]=  

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.
  
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
f[1], f[2], ....
By default, the new parameters generated by Reduce are named C[1], C[2], ....
Out[49]=  

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
  
"BranchLinearDiophantine"  True  whether Reduce should use a branchandbound 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 
"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  the maximal number of points at which the modular sieve method evaluates the system 
ReduceOptions group options affecting the behavior of Reduce, Resolve, and FindInstance for Diophantine polynomial systems.
BranchLinearDiophantine
The value of the system option
BranchLinearDiophantine specifies which variant of the algorithm should be used in
the final stage of solving bounded linear systems. Neither variant seems to be clearly better. For some examples the hybrid method combining a branchandbound type algorithm and a simple recursive enumeration is faster; for other examples the simple recursive enumeration alone is faster. The hybrid method seems to be more robust for badly conditioned problems, hence it is the default method.
This finds integer points in a long, narrow fourdimensional simplex using the default hybrid method.
Out[52]=  

This sets the value of the system option BranchLinearDiophantine to False. 
Here the simple recursive enumeration method is used, and for this badly conditioned problem it is several times slower.
Out[54]=  

This resets the value of the system option BranchLinearDiophantine to the default value. 
Here are solutions of a system of two randomly generated equations eqns and three randomly generated inequalities ineqs in seven variables inside a simplex bounded by bds.
Out[64]=  

For this system the nondefault simple recursion method is faster.
Out[66]=  

Here is a random system very similar to the previous one, except that it contains one more variable and the righthand side of the last of bds is changed from 100 to 200. However, for this system the default method is faster.
Out[76]=  

The nondefault method is slower for this system.
Out[78]=  

This resets the value of the system option BranchLinearDiophantine to the default value. 
DiscreteSolutionBound
The value of the system option
DiscreteSolutionBound specifies whether integer solutions in a real interval
a≤x≤b should be enumerated explicitly or represented implicitly as
x a≤x≤b. 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 0≤x<10. Reduce writes them out explicitly.
Out[80]=  

There are 11 integers in the real interval 0≤x<1001^{1/3}. Reduce represents them implicitly.
Out[81]=  

This increases the DiscreteSolutionBound to 11. 
This resets DiscreteSolutionBound to the default value. 
The value of
DiscreteSolutionBound also affects
the solving of bounded linear systems.
ExhaustiveSearchMaxPoints
The system option
ExhaustiveSearchMaxPoints 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
{1000, 10000}). 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.
With the default setting of ExhaustiveSearchMaxPoints, Reduce is unable to solve this equation.
Out[85]=  

This increases the value of the second element of ExhaustiveSearchMaxPoints to 100000. 
With the default setting of ExhaustiveSearchMaxPoints, Reduce solves this equation using the method for solving Pell equations.
Out[88]=  

Increasing the first element of ExhaustiveSearchMaxPoints to 10^{6} makes Reduce use the exhaustive search first. In this example the search is much slower than the Pell equation solver.
Out[90]=  

For this equation the Pell equation solver is slower than the exhaustive search.
Out[92]=  

The exhaustive search is faster here.
Out[93]=  

This resets ExhaustiveSearchMaxPoints to the default value. 
LatticeReduceDiophantine
The value of the system option
LatticeReduceDiophantine 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.
This finds the only two integer points in a triangle whose projections on both axes have sizes greater than a but whose projection on the line x+5000y=0 has size one.
Out[96]=  

This sets the value of the system option LatticeReduceDiophantine to False. 
The nondefault method is much slower for this system, and the speed difference grows with a.
Out[98]=  

Here is a system that contains a set of simple inequalities bds, which bound solutions to a reasonably small size polyhedron, combined with a set of relatively complicated inequalities ineqs. For such systems, using LatticeReduce tends to increase the timing.
Out[105]=  

The nondefault method is faster for this system.
Out[107]=  

This resets LatticeReduceDiophantine to the default value. 
MaxFrobeniusGraph
The system option
MaxFrobeniusGraph 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.
To find a solution of a Frobenius equation with the smallest coefficient larger than 10^{6}, FindInstance by default uses the general method for solving bounded linear systems. For this example the method is relatively slow but uses little memory. The kernel has been restarted to show the memory usage by the current example.
Out[4]=  
Out[5]=  

This increases the value of MaxFrobeniusGraph to 10^{7}. 
Now FindInstance uses the method based on the computation of the Frobenius graph. It finds the solution faster, but it uses more memory.
Out[7]=  
Out[8]=  

This resets MaxFrobeniusGraph to the default value. 
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 10,000.
With the default setting of SieveMaxPoints, FindInstance is unable to find a solution for this equation.
Out[10]=  

Increasing the number of SieveMaxPoints to one million allows FindInstance to find a solution.
Out[12]=  

This resets SieveMaxPoints to the default value. 
References
[1] C. Goldbach, Letter to L. Euler, June 7, 1742. http://www.mathstat.dal.ca/~joerg/pic/gletter.jpg.
[2] Matiyasevich Yu. V. "The Diophantineness of Enumerable Sets (Russian)" Dokl. Akad. Nauk SSSR 191 (1970): 279282. English translation: Soviet Math. Dokl. 11 (1970): 354358
[3] Contejean E. and H. Devie. "An Efficient Incremental Algorithm for Solving Systems of Linear Diophantine Equations" Information and Computation 113 (1994): 143172
[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): 2130
[5] Strzebonski A. "An Improved Algorithm for Diophantine Equations in One Variable" Paper presented at the ACA 2002 Session on SymbolicNumerical 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 n=f u^{2}+g v^{2} in Coprime Integers u and v" Mathematics of Computation 55 (1990): 327343
[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): 239256