Navigation

The word “Diophantine” comes with the name Diophantus, a mathematician lived
in the great city of Alexandria sometime around 250 AD. Often referred to as
the “father of Algebra”, Diophantus in his famous work “Arithmetica”
presented 150 problems that marked the early beginnings of number theory, the
field of study about integers and their properties. Diophantine equations play
a central and an important part in number theory.

We call a “Diophantine equation” to an equation of the form,
\(f(x_1, x_2, \ldots x_n) = 0\) where \(n \geq 2\) and \(x_1, x_2, \ldots x_n\) are
integer variables. If we can find \(n\) integers \(a_1, a_2, \ldots a_n\) such that
\(x_1 = a_1, x_2 = a_2, \ldots x_n = a_n\) satisfies the above equation, we say
that the equation is solvable. You can read more about Diophantine equations in
[1] and [2].

Currently, following five types of Diophantine equations can be solved using
diophantine() and other helper functions of
the Diophantine module.

When an equation is given to diophantine(),
it factors the equation(if possible) and solves the equation given by each
factor by calling diop_solve() separately.
Then all the results are combined using merge_solution().

Before we start solving the equations, we need to define the variables.

>>> fromsympyimportsymbols>>> x,y,z=symbols("x, y, z",integer=True)

Let’s start by solving the easiest type of Diophantine equations, i.e. linear
Diophantine equations. Let’s solve \(2x + 3y = 5\). Note that although we
write the equation in the above form, when we input the equation to any of the
functions in Diophantine module, it needs to be in the form \(eq = 0\).

>>> diophantine(2*x+3*y-5)set([(3*t - 5, -2*t + 5)])

Note that stepping one more level below the highest API, we can solve the very
same equation by calling diop_solve().

Note that except for the highest level API, in case of no solutions, a tuple of
\(None\) are returned. Size of the tuple is the same as the number of variables.
Also, one can specifically set the parameter to be used in the solutions by
passing a customized parameter. Consider the following example.

Please note that for the moment, user can set the parameter only for linear
Diophantine equations and binary quadratic equations.

Let’s try solving a binary quadratic equation which is an equation with two
variables and has a degree of two. Before trying to solve these equations, an
idea about various cases associated with the equation would help a lot. Please
refer [3] and [4] for detailed analysis of different cases and the nature
of the solutions. Let us define \(\Delta = b^2 - 4ac\) w.r.t. the binary quadratic
\(ax^2 + bxy + cy^2 + dx + ey + f = 0\).

When \(\Delta < 0\), there are either no solutions or only a finite number of solutions.

>>> diophantine(x**2-4*x*y+8*y**2-3*x+7*y-5)set([(2, 1), (5, 1)])

In the above equation \(\Delta = (-4)^2 - 4*1*8 = -16\) and hence only a finite
number of solutions exist.

When \(\Delta = 0\) we might have either no solutions or parameterized solutions.

The most interesting case is when \(\Delta > 0\) and it is not a perfect square.
In this case, the equation has either no solutions or an infinte number of
solutions. Consider the below cases where \(\Delta = 8\).

Here \(n\) is an integer. Although x_n and y_n may not look like
integers, substituting in specific values for n (and simplifying) shows that they
are. For example consider the following example where we set n equal to 9.

gives the equation \(X^2 -5Y^2 = 920\). Values of \(A\) and \(B\) are as belows.

>>> AMatrix([[1/10, 3/10],[ 0, 1/5]])>>> BMatrix([[ 1/5],[-11/5]])

We can solve an equation of the form \(X^2 - DY^2 = N\) by passing \(D\) and \(N\) to
diop_DN()

>>> diop_DN(5,920)[]

Unfortunately, our equation does not have solutions.

Now let’s turn to homogeneous ternary quadratic equations. These equations are
of the form \(ax^2 + by^2 + cz^2 + dxy + eyz + fzx = 0\). These type of equations
either have infinitely many solutions or no solutions (except the obvious
solution (0, 0, 0))

diop_ternary_quadratic() first converts the
given equation to an equivalent equation of the form \(w^2 = AX^2 + BY^2\) and
then it uses descent() to solve the latter
equation. You can refer to the docs of
transformation_to_normal() to find more on
this. The equation \(w^2 = AX^2 + BY^2\) can be solved more easily by using the
Aforementioned descent().

Simplify the solution procedure of diophantine equation eq by
converting it into a product of terms which should equal zero.

For example, when solving, \(x^2 - y^2 = 0\) this is treated as
\((x + y)(x - y) = 0\) and \(x+y = 0\) and \(x-y = 0\) are solved independently
and combined. Each term is solved by calling diop_solve().

Output of diophantine() is a set of tuples. Each tuple represents a
solution of the input equation. In a tuple, solution for each variable is
listed according to the alphabetic order of input variables. i.e. if we have
an equation with two variables \(a\) and \(b\), first element of the tuple will
give the solution for \(a\) and the second element will give the solution for
\(b\).

Returns a tuple containing the type of the diophantine equation along with
the variables(free symbols) and their coefficients. Variables are returned
as a list and coefficients are returned as a dict with the key being the
respective term and the constant term is keyed to Integer(1). Type is an
element in the set {“linear”, “binary_quadratic”, “general_pythagorean”,
“homogeneous_ternary_quadratic”, “univariate”, “general_sum_of_squares”}

i.e. equations of the form \(Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0\). Returns a
set containing the tuples \((x, y)\) which contains the solutions. If there
are no solutions then \((None, None)\) is returned.

Mainly concerned in the case \(D > 0, D\) is not a perfect square, which is
the same as generalized Pell equation. To solve the generalized Pell
equation this function Uses LMM algorithm. Refer [R357] for more details on
the algorithm.
Returns one solution for each class of the solutions. Other solutions of
the class can be constructed according to the values of D and N.
Returns a list containing the solution tuples \((x, y)\).

The output can be interpreted as follows: There are three fundamental
solutions to the equation \(x^2 - 13y^2 = -4\) given by (3, 1), (393, 109)
and (36, 10). Each tuple is in the form (x, y), i. e solution (3, 1) means
that \(x = 3\) and \(y = 1\).

Uses the algorithm due to Cornacchia. The method only finds primitive
solutions, i.e. ones with \(\gcd(x, y) = 1\). So this method can’t be used to
find the solutions of \(x^2 + y^2 = 20\) since the only solution to former is
\((x,y) = (4, 2)\) and it is not primitive. When ` a = b = 1`, only the
solutions with \(x \geq y\) are found. For more details, see the References.

Mainly concerned with the generalized Pell equation which is the case when
\(D > 0, D\) is not a perfect square. For more information on the case refer
[R360]. Let \((t, u)\) be the minimal positive solution of the equation
\(x^2 - Dy^2 = 1\). Then this method requires
\(\sqrt{\frac{\mid N \mid (t \pm 1)}{2D}}\) to be small.

This is used to solve the general quadratic equation by transforming it to
the latter form. Refer [R361] for more detailed information on the
transformation. This function returns a tuple (A, B) where A is a 2 X 2
matrix and B is a 2 X 1 matrix such that,

A, B returned are such that Transpose((x y)) = A * Transpose((X Y)) + B.
Substituting these values for \(x\) and \(y\) and a bit of simplifying work
will give an equation of the form \(x^2 - Dy^2 = N\).

general_sum_of_squares(eq,limit) : Here eq is an expression which
is assumed to be zero. Also, eq should be in the form,
\(x_{1}^2 + x_{2}^2 + . . . + x_{n}^2 - k = 0\). At most limit number of
solutions are returned.

Details

When \(n = 3\) if \(k = 4^a(8m + 7)\) for some \(a, m \in Z\) then there will be
no solutions. Refer [R364] for more details.

Returns a generator that can be used to generate partitions of an integer
\(n\).

A partition of \(n\) is a set of positive integers which add upto \(n\). For
example, partitions of 3 are 3 , 1 + 2, 1 + 1+ 1. A partition is returned
as a tuple. If k equals None, then all possible partitions are returned
irrespective of their size, otherwise only the partitions of size k are
returned. If there are no partions of \(n\) with size \(k\) then an empty tuple
is returned. If the zero parameter is set to True then a suitable
number of zeros are added at the end of every partition of size less than
k.

zero parameter is considered only if k is not None. When the
partitions are over, the last \(next()\) call throws the StopIteration
exception, so this function should always be used inside a try - except
block.

This is used to construct the full solution from the solutions of sub
equations.

For example when solving the equation \((x - y)(x^2 + y^2 - z^2) = 0\),
solutions for each of the equations \(x-y = 0\) and \(x^2 + y^2 - z^2\) are
found independently. Solutions for \(x - y = 0\) are \((x, y) = (t, t)\). But
we should introduce a value for z when we output the solution for the
original equation. This function converts \((t, t)\) into \((t, t, n_{1})\)
where \(n_{1}\) is an integer parameter.

There are six sequences of integers defined related to the continued
fraction representation of \(\frac{P + \sqrt{D}}{Q}\), namely {\(P_{i}\)},
{\(Q_{i}\)}, {\(a_{i}\)},{\(A_{i}\)}, {\(B_{i}\)}, {\(G_{i}\)}. PQa() Returns
these values as a 6-tuple in the same order as mentioned above. Refer [R368]
for more detailed information.

PQa(P_0,Q_0,D): P_0, Q_0 and D are integers corresponding
to \(P_{0}\), \(Q_{0}\) and \(D\) in the continued fraction
\(\frac{P_{0} + \sqrt{D}}{Q_{0}}\).
Also it’s assumed that \(P_{0}^2 == D mod(|Q_{0}|)\) and \(D\) is square free.

Returns True if two solutions \((u, v)\) and \((r, s)\) of \(x^2 - Dy^2 = N\)
belongs to the same equivalence class and False otherwise.

Two solutions \((u, v)\) and \((r, s)\) to the above equation fall to the same
equivalence class iff both \((ur - Dvs)\) and \((us - vr)\) are divisible by
\(N\). See reference [R369]. No test is performed to test whether \((u, v)\) and
\((r, s)\) are actually solutions to the equation. User should take care of
this.

Here the coefficients \(a\), \(b\), and \(c\) should be non zero. Otherwise the
equation will be a quadratic binary or univariate equation. If solvable,
returns a tuple \((x, y, z)\) that satisifes the given equation. If the
equation does not have integer solutions, \((None, None, None)\) is returned.