Diophantine

Diophantine equations

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.

  • Linear Diophantine equations: \(a_1x_1 + a_2x_2 + \ldots + a_nx_n = b\).
  • General binary quadratic equation: \(ax^2 + bxy + cy^2 + dx + ey + f = 0\)
  • Homogeneous ternary quadratic equation: \(ax^2 + by^2 + cz^2 + dxy + eyz + fzx = 0\)
  • Extended Pythagorean equation: \(a_{1}x_{1}^2 + a_{2}x_{2}^2 + \ldots + a_{n}x_{n}^2 = a_{n+1}x_{n+1}^2\)
  • General sum of squares: \(x_{1}^2 + x_{2}^2 + \ldots + x_{n}^2 = k\)

Module structure

This module contains diophantine() and helper functions that are needed to solve certain Diophantine equations. It’s structured in the following manner.

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

diop_solve() internally uses classify_diop() to find the type of the equation(and some other details) given to it and then calls the appropriate solver function based on the type returned. For example, if classify_diop() returned “linear” as the type of the equation, then diop_solve() calls diop_linear() to solve the equation.

Each of the functions, diop_linear(), diop_quadratic(), diop_ternary_quadratic(), diop_general_pythagorean() and diop_general_sum_of_squares() solves a specific type of equations and the type can be easily guessed by it’s name.

Apart from these functions, there are a considerable number of other functions in the “Diophantine Module” and all of them are listed under User functions and Internal functions.

Tutorial

First, let’s import the highest API of the Diophantine module.

>>> from diofant.solvers.diophantine import diophantine

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

>>> from diofant import symbols
>>> x, y, z, t, p, q = symbols("x, y, z, t, p, q", integer=True)
>>> t1, t2, t3, t4, t5 = symbols("t1:6", 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)
{(3*t_0 - 5, -2*t_0 + 5)}

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

>>> from diofant.solvers.diophantine import diop_solve
>>> diop_solve(2*x + 3*y - 5)
(3*t_0 - 5, -2*t_0 + 5)

Note that it returns a tuple rather than a set. diophantine() always return a set of tuples. But diop_solve() may return a single tuple or a set of tuples depending on the type of the equation given.

We can also solve this equation by calling diop_linear(), which is what diop_solve() calls internally.

>>> from diofant.solvers.diophantine import diop_linear
>>> diop_linear(2*x + 3*y - 5)
(3*t_0 - 5, -2*t_0 + 5)

If the given equation has no solutions then the outputs will look like below.

>>> diophantine(2*x + 4*y - 3)
set()
>>> diop_solve(2*x + 4*y - 3)
(None, None)
>>> diop_linear(2*x + 4*y - 3)
(None, None)

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:

>>> m = symbols("m", integer=True)
>>> diop_solve(2*x + 3*y - 5, m)
(3*m_0 - 5, -2*m_0 + 5)

For linear Diophantine equations, the customized parameter is the prefix used for each free variable in the solution. Consider the following example:

>>> diop_solve(2*x + 3*y - 5*z + 7, m)
(m_0, -9*m_0 - 5*m_1 - 14, -5*m_0 - 3*m_1 - 7)

In the solution above, m_0 and m_1 are independent free variables.

Please note that for the moment, users 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)
{(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.

>>> diophantine(3*x**2 - 6*x*y + 3*y**2 - 3*x + 7*y - 5)
set()
>>> diophantine(x**2 - 4*x*y + 4*y**2 - 3*x + 7*y - 5)
{(-2*t**2 - 7*t + 10, -t**2 - 3*t + 5)}
>>> diophantine(x**2 + 2*x*y + y**2 - 3*x - 3*y)
{(t_0, -t_0), (t_0, -t_0 + 3)}

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\).

>>> diophantine(x**2 - 4*x*y + 2*y**2 - 3*x + 7*y - 5)
set()
>>> from diofant import sqrt
>>> n = symbols("n", integer=True)
>>> s = diophantine(x**2 -  2*y**2 - 2*x - 4*y, n)
>>> x_1, y_1 = s.pop()
>>> x_2, y_2 = s.pop()
>>> x_n = -(-2*sqrt(2) + 3)**n/2 + sqrt(2)*(-2*sqrt(2) + 3)**n/2 - sqrt(2)*(2*sqrt(2) + 3)**n/2 - (2*sqrt(2) + 3)**n/2 + 1
>>> x_1 == x_n or x_2 == x_n
True
>>> y_n = -sqrt(2)*(-2*sqrt(2) + 3)**n/4 + (-2*sqrt(2) + 3)**n/2 + sqrt(2)*(2*sqrt(2) + 3)**n/4 + (2*sqrt(2) + 3)**n/2 - 1
>>> y_1 == y_n or y_2 == y_n
True

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.

>>> from diofant import simplify
>>> simplify(x_n.subs({n: 9}))
-9369318

Any binary quadratic of the form \(ax^2 + bxy + cy^2 + dx + ey + f = 0\) can be transformed to an equivalent form \(X^2 - DY^2 = N\).

>>> from diofant.solvers.diophantine import find_DN, diop_DN, transformation_to_DN
>>> find_DN(x**2 - 3*x*y + y**2 - 7*x + 5*y - 3)
(5, 920)

So, the above equation is equivalent to the equation \(X^2 - 5Y^2 = 920\) after a linear transformation. If we want to find the linear transformation, we can use transformation_to_DN()

>>> A, B = transformation_to_DN(x**2 - 3*x*y + y**2 - 7*x + 5*y - 3)

Here A is a 2 X 2 matrix and B is a 2 X 1 matrix such that the transformation

\[\begin{split}\begin{bmatrix} X\\Y \end{bmatrix} = A \begin{bmatrix} x\\y \end{bmatrix} + B\end{split}\]

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

>>> A
Matrix([
[1/10, 3/10],
[   0,  1/5]])
>>> B
Matrix([
[  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))

>>> diophantine(3*x**2 + 4*y**2 - 5*z**2 + 4*x*y + 6*y*z + 7*z*x)
set()
>>> diophantine(3*x**2 + 4*y**2 - 5*z**2 + 4*x*y - 7*y*z + 7*z*x)
{(-16*p**2 + 28*p*q + 20*q**2, 3*p**2 + 38*p*q - 25*q**2, 4*p**2 - 24*p*q + 68*q**2)}

If you are only interested about a base solution rather than the parameterized general solution (to be more precise, one of the general solutions), you can use diop_ternary_quadratic().

>>> from diofant.solvers.diophantine import diop_ternary_quadratic
>>> diop_ternary_quadratic(3*x**2 + 4*y**2 - 5*z**2 + 4*x*y - 7*y*z + 7*z*x)
(-4, 5, 1)

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

>>> from diofant.solvers.diophantine import descent
>>> descent(3, 1) # solves the equation w**2 = 3*Y**2 + Z**2
(1, 0, 1)

Here the solution tuple is in the order (w, Y, Z)

The extended Pythagorean equation, \(a_{1}x_{1}^2 + a_{2}x_{2}^2 + \ldots + a_{n}x_{n}^2 = a_{n+1}x_{n+1}^2\) and the general sum of squares equation, \(x_{1}^2 + x_{2}^2 + \ldots + x_{n}^2 = k\) can also be solved using the Diophantine module.

>>> from diofant.abc import a, b, c, d, e, f
>>> diophantine(9*a**2 + 16*b**2 + c**2 + 49*d**2 + 4*e**2 - 25*f**2)
{(70*t1**2 + 70*t2**2 + 70*t3**2 + 70*t4**2 - 70*t5**2, 105*t1*t5, 420*t2*t5, 60*t3*t5, 210*t4*t5, 42*t1**2 + 42*t2**2 + 42*t3**2 + 42*t4**2 + 42*t5**2)}

function diop_general_pythagorean() can also be called directly to solve the same equation. This is true about the general sum of squares too. Either you can call diop_general_pythagorean() or use the high level API.

>>> diophantine(a**2 + b**2 + c**2 + d**2 + e**2 + f**2 - 112)
{(8, 4, 4, 4, 0, 0)}

If you want to get a more thorough idea about the the Diophantine module please refer to the following blog.

http://thilinaatdiofant.wordpress.com/

References

[1]Andreescu, Titu. Andrica, Dorin. Cucurezeanu, Ion. An Introduction to Diophantine Equations
[2]Diophantine Equation, Wolfram Mathworld, [online]. Available: http://mathworld.wolfram.com/DiophantineEquation.html
[3]Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0,[online], Available: http://www.alpertron.com.ar/METHODS.HTM
[4]Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online], Available: http://www.jpr2718.org/ax2p.pdf

User Functions

These are functions that are imported into the global namespace with from diofant import *. These functions are intended for use by ordinary users of Diofant.

diophantine

diofant.solvers.diophantine.diophantine(eq, param=Symbol('t', integer=True))[source]

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\).

Parameters:

eq : Relational or Expr

an equation (to be solved)

t : Symbol, optional

the parameter to be used in the solution.

Examples

>>> from diofant.solvers.diophantine import diophantine
>>> from diofant.abc import x, y, z
>>> diophantine(x**2 - y**2)
{(-t_0, -t_0), (t_0, -t_0)}

diop_solve

diofant.solvers.diophantine.diop_solve(eq, param=Symbol('t', integer=True))[source]

Solves the diophantine equation eq.

Similar to diophantine() but doesn’t try to factor eq as latter does. Uses classify_diop() to determine the type of the eqaution and calls the appropriate solver function.

Parameters:

eq : Expr

an expression, which is assumed to be zero.

t : Symbol, optional

a parameter, to be used in the solution.

Examples

>>> from diofant.solvers.diophantine import diop_solve
>>> from diofant.abc import x, y, z, w
>>> diop_solve(2*x + 3*y - 5)
(3*t_0 - 5, -2*t_0 + 5)
>>> diop_solve(4*x + 3*y -4*z + 5)
(t_0, -4*t_1 + 5, t_0 - 3*t_1 + 5)
>>> diop_solve(x + 3*y - 4*z + w -6)
(t_0, t_0 + t_1, -2*t_0 - 3*t_1 - 4*t_2 - 6, -t_0 - 2*t_1 - 3*t_2 - 6)
>>> diop_solve(x**2 + y**2 - 5)
{(-2, -1), (-2, 1), (2, -1), (2, 1)}

classify_diop

diofant.solvers.diophantine.classify_diop(eq)[source]

Helper routine used by diop_solve() to find the type of the eq etc.

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”}

Parameters:

eq : Expr

an expression, which is assumed to be zero.

Examples

>>> from diofant.solvers.diophantine import classify_diop
>>> from diofant.abc import x, y, z, w, t
>>> classify_diop(4*x + 6*y - 4)
([x, y], {1: -4, x: 4, y: 6}, 'linear')
>>> classify_diop(x + 3*y -4*z + 5)
([x, y, z], {1: 5, x: 1, y: 3, z: -4}, 'linear')
>>> classify_diop(x**2 + y**2 - x*y + x + 5)
([x, y], {1: 5, x: 1, x**2: 1, y: 0, y**2: 1, x*y: -1}, 'binary_quadratic')

diop_linear

diofant.solvers.diophantine.diop_linear(eq, param=Symbol('t', integer=True))[source]

Solves linear diophantine equations.

A linear diophantine equation is an equation of the form \(a_{1}x_{1} + a_{2}x_{2} + .. + a_{n}x_{n} = 0\) where \(a_{1}, a_{2}, ..a_{n}\) are integer constants and \(x_{1}, x_{2}, ..x_{n}\) are integer variables.

Parameters:

eq : Expr

is a linear diophantine equation which is assumed to be zero.

param : Symbol, optional

is the parameter to be used in the solution.

Examples

>>> from diofant.solvers.diophantine import diop_linear
>>> from diofant.abc import x, y, z
>>> from diofant import Integer
>>> diop_linear(2*x - 3*y - 5)  # solves equation 2*x - 3*y -5 = 0
(-3*t_0 - 5, -2*t_0 - 5)

Here x = -3*t_0 - 5 and y = -2*t_0 - 5

>>> diop_linear(2*x - 3*y - 4*z -3)
(t_0, -6*t_0 - 4*t_1 + 3, 5*t_0 + 3*t_1 - 3)

base_solution_linear

diofant.solvers.diophantine.base_solution_linear(c, a, b, t=None)[source]

Return the base solution for a linear diophantine equation with two variables.

Used by diop_linear() to find the base solution of a linear Diophantine equation. If t is given then the parametrized solution is returned.

base_solution_linear(c, a, b, t): a, b, c are coefficients in \(ax + by = c\) and t is the parameter to be used in the solution.

Examples

>>> from diofant.solvers.diophantine import base_solution_linear
>>> from diofant.abc import t
>>> base_solution_linear(5, 2, 3)  # equation 2*x + 3*y = 5
(-5, 5)
>>> base_solution_linear(0, 5, 7)  # equation 5*x + 7*y = 0
(0, 0)
>>> base_solution_linear(5, 2, 3, t)  # equation 2*x + 3*y = 5
(3*t - 5, -2*t + 5)
>>> base_solution_linear(0, 5, 7, t)  # equation 5*x + 7*y = 0
(7*t, -5*t)

diop_quadratic

diofant.solvers.diophantine.diop_quadratic(eq, param=Symbol('t', integer=True))[source]

Solves quadratic diophantine equations.

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.

Parameters:

eq : Expr

should be a quadratic bivariate expression which is assumed to be zero.

param : Symbol, optional

is a parameter to be used in the solution.

References

[R471]Methods to solve Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0,[online], Available: http://www.alpertron.com.ar/METHODS.HTM
[R472]Solving the equation ax^2+ bxy + cy^2 + dx + ey + f= 0, [online], Available: http://www.jpr2718.org/ax2p.pdf

Examples

>>> from diofant.abc import x, y, t
>>> from diofant.solvers.diophantine import diop_quadratic
>>> diop_quadratic(x**2 + y**2 + 2*x + 2*y + 2, t)
{(-1, -1)}

diop_DN

diofant.solvers.diophantine.diop_DN(D, N, t=Symbol('t', integer=True))[source]

Solves the equation \(x^2 - Dy^2 = N\).

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 [R473] 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)\).

Parameters:

D, N : Integer

correspond to D and N in the equation.

t : Symbol, optional

is the parameter to be used in the solutions.

References

[R473](1, 2) Solving the generalized Pell equation x**2 - D*y**2 = N, John P. Robertson, July 31, 2004, Pages 16 - 17. [online], Available: http://www.jpr2718.org/pell.pdf

Examples

>>> from diofant.solvers.diophantine import diop_DN
>>> diop_DN(13, -4)  # Solves equation x**2 - 13*y**2 = -4
[(3, 1), (393, 109), (36, 10)]

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\).

>>> diop_DN(986, 1)  # Solves equation x**2 - 986*y**2 = 1
[(49299, 1570)]

cornacchia

diofant.solvers.diophantine.cornacchia(a, b, m)[source]

Solves \(ax^2 + by^2 = m\) where \(\gcd(a, b) = 1 = gcd(a, m)\) and \(a, b > 0\).

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.

References

[R474]
  1. Nitaj, “L’algorithme de Cornacchia”
[R475]Solving the diophantine equation ax**2 + by**2 = m by Cornacchia’s method, [online], Available: http://www.numbertheory.org/php/cornacchia.html

Examples

>>> from diofant.solvers.diophantine import cornacchia
>>> cornacchia(2, 3, 35)  # equation 2x**2 + 3y**2 = 35
{(2, 3), (4, 1)}
>>> cornacchia(1, 1, 25)  # equation x**2 + y**2 = 25
{(4, 3)}

diop_bf_DN

diofant.solvers.diophantine.diop_bf_DN(D, N, t=Symbol('t', integer=True))[source]

Uses brute force to solve the equation, \(x^2 - Dy^2 = N\).

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 [R476]. 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.

Parameters:

D, N : Integer

correspond to D and N in the equation.

t : Symbol, optional

is the parameter to be used in the solutions.

References

[R476](1, 2) Solving the generalized Pell equation x**2 - D*y**2 = N, John P. Robertson, July 31, 2004, Page 15. http://www.jpr2718.org/pell.pdf

Examples

>>> from diofant.solvers.diophantine import diop_bf_DN
>>> diop_bf_DN(13, -4)
[(3, 1), (-3, 1), (36, 10)]
>>> diop_bf_DN(986, 1)
[(49299, 1570)]

transformation_to_DN

diofant.solvers.diophantine.transformation_to_DN(eq)[source]

This function transforms general quadratic, \(ax^2 + bxy + cy^2 + dx + ey + f = 0\) to more easy to deal with \(X^2 - DY^2 = N\) form.

This is used to solve the general quadratic equation by transforming it to the latter form. Refer [R477] 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,

Transpose([x y]) = A * Transpose([X Y]) + B

Parameters:

eq : Expr

the quadratic expression to be transformed.

References

[R477](1, 2) Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0, John P.Robertson, May 8, 2003, Page 7 - 11. http://www.jpr2718.org/ax2p.pdf

Examples

>>> from diofant.abc import x, y
>>> from diofant.solvers.diophantine import transformation_to_DN
>>> from diofant.solvers.diophantine import classify_diop
>>> A, B = transformation_to_DN(x**2 - 3*x*y - y**2 - 2*y + 1)
>>> A
Matrix([
[1/26, 3/26],
[   0, 1/13]])
>>> B
Matrix([
[-6/13],
[-4/13]])

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\).

>>> from diofant.abc import X, Y
>>> from diofant import Matrix, simplify, Subs
>>> u = (A*Matrix([X, Y]) + B)[0]  # Transformation for x
>>> u
X/26 + 3*Y/26 - 6/13
>>> v = (A*Matrix([X, Y]) + B)[1]  # Transformation for y
>>> v
Y/13 - 4/13

Next we will substitute these formulas for \(x\) and \(y\) and do simplify().

>>> eq = simplify(Subs(x**2 - 3*x*y - y**2 - 2*y + 1, (x, y), (u, v)).doit())
>>> eq
X**2/676 - Y**2/52 + 17/13

By multiplying the denominator appropriately, we can get a Pell equation in the standard form.

>>> eq * 676
X**2 - 13*Y**2 + 884

If only the final equation is needed, find_DN() can be used.

find_DN

diofant.solvers.diophantine.find_DN(eq)[source]

This function returns a tuple, \((D, N)\) of the simplified form, \(x^2 - Dy^2 = N\), corresponding to the general quadratic, \(ax^2 + bxy + cy^2 + dx + ey + f = 0\).

Solving the general quadratic is then equivalent to solving the equation \(X^2 - DY^2 = N\) and transforming the solutions by using the transformation matrices returned by transformation_to_DN().

Parameters:

eq : Expr

is the quadratic expression to be transformed.

References

[R478]Solving the equation ax^2 + bxy + cy^2 + dx + ey + f = 0, John P.Robertson, May 8, 2003, Page 7 - 11. http://www.jpr2718.org/ax2p.pdf

Examples

>>> from diofant.abc import x, y
>>> from diofant.solvers.diophantine import find_DN
>>> find_DN(x**2 - 3*x*y - y**2 - 2*y + 1)
(13, -884)

Interpretation of the output is that we get \(X^2 -13Y^2 = -884\) after transforming \(x^2 - 3xy - y^2 - 2y + 1\) using the transformation returned by transformation_to_DN().

diop_ternary_quadratic

diofant.solvers.diophantine.diop_ternary_quadratic(eq)[source]

Solves the general quadratic ternary form, \(ax^2 + by^2 + cz^2 + fxy + gyz + hxz = 0\).

Parameters:

eq : Expr

should be an homogeneous expression of degree two in three variables and it is assumed to be zero.

Returns:

tuple

which is a base solution for the above equation. If there are no solutions, (None, None, None) is returned.

Examples

>>> from diofant.abc import x, y, z
>>> from diofant.solvers.diophantine import diop_ternary_quadratic
>>> diop_ternary_quadratic(x**2 + 3*y**2 - z**2)
(1, 0, 1)
>>> diop_ternary_quadratic(4*x**2 + 5*y**2 - z**2)
(1, 0, 2)
>>> diop_ternary_quadratic(45*x**2 - 7*y**2 - 8*x*y - z**2)
(28, 45, 105)
>>> diop_ternary_quadratic(x**2 - 49*y**2 - z**2 + 13*z*y -8*x*y)
(9, 1, 5)

square_factor

diofant.solvers.diophantine.square_factor(a)[source]

Returns an integer \(c\) s.t. \(a = c^2k, \ c,k \in Z\). Here \(k\) is square free.

Examples

>>> from diofant.solvers.diophantine import square_factor
>>> square_factor(24)
2
>>> square_factor(36)
6
>>> square_factor(1)
1

descent

diofant.solvers.diophantine.descent(A, B)[source]

Lagrange’s \(descent()\) with lattice-reduction to find solutions to \(x^2 = Ay^2 + Bz^2\).

Here \(A\) and \(B\) should be square free and pairwise prime. Always should be called with suitable A and B so that the above equation has solutions.

This is more faster than the normal Lagrange’s descent algorithm because the gaussian reduction is used.

References

[R479]Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin, Mathematics of Computation, Volume 00, Number 0.

Examples

>>> from diofant.solvers.diophantine import descent
>>> descent(3, 1)  # x**2 = 3*y**2 + z**2
(1, 0, 1)

\((x, y, z) = (1, 0, 1)\) is a solution to the above equation.

>>> descent(41, -113)
(-16, -3, 1)

diop_general_pythagorean

diofant.solvers.diophantine.diop_general_pythagorean(eq, param=Symbol('m', integer=True))[source]

Solves the general pythagorean equation, \(a_{1}^2x_{1}^2 + a_{2}^2x_{2}^2 + . . . + a_{n}^2x_{n}^2 - a_{n + 1}^2x_{n + 1}^2 = 0\).

Returns a tuple which contains a parametrized solution to the equation, sorted in the same order as the input variables.

Parameters:

eq : Expr

is a general pythagorean equation which is assumed to be zero

param : Symbol, optional

is the base parameter used to construct other parameters by subscripting.

Examples

>>> from diofant.solvers.diophantine import diop_general_pythagorean
>>> from diofant.abc import a, b, c, d, e
>>> diop_general_pythagorean(a**2 + b**2 + c**2 - d**2)
(m1**2 + m2**2 - m3**2, 2*m1*m3, 2*m2*m3, m1**2 + m2**2 + m3**2)
>>> diop_general_pythagorean(9*a**2 - 4*b**2 + 16*c**2 + 25*d**2 + e**2)
(10*m1**2  + 10*m2**2  + 10*m3**2 - 10*m4**2, 15*m1**2  + 15*m2**2  + 15*m3**2  + 15*m4**2, 15*m1*m4, 12*m2*m4, 60*m3*m4)

diop_general_sum_of_squares

diofant.solvers.diophantine.diop_general_sum_of_squares(eq, limit=1)[source]

Solves the equation \(x_{1}^2 + x_{2}^2 + . . . + x_{n}^2 - k = 0\).

Returns at most limit number of solutions. Currently there is no way to set limit using higher level API’s like diophantine() or diop_solve() but that will be fixed soon.

Parameters:

eq : Expr

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\).

limit : int, optional

upper limit (the default is 1) for number of solutions returned.

Notes

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

References

[R480](1, 2) Representing an Integer as a sum of three squares, [online], Available: http://www.proofwiki.org/wiki/Integer_as_Sum_of_Three_Squares

Examples

>>> from diofant.solvers.diophantine import diop_general_sum_of_squares
>>> from diofant.abc import a, b, c, d, e, f
>>> diop_general_sum_of_squares(a**2 + b**2 + c**2 + d**2 + e**2 - 2345) == {(0, 48, 5, 4, 0)}
True

partition

diofant.solvers.diophantine.partition(n, k=None, zeros=False)[source]

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

Parameters:

n : int

is a positive integer

k : int, optional

is the size of the partition which is also positive integer. The default is None.

zeros : boolean, optional

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.

References

[R481]Generating Integer Partitions, [online], Available: http://jeromekelleher.net/partitions.php

Examples

>>> from diofant.solvers.diophantine import partition
>>> f = partition(5)
>>> next(f)
(1, 1, 1, 1, 1)
>>> next(f)
(1, 1, 1, 2)
>>> g = partition(5, 3)
>>> next(g)
(3, 1, 1)
>>> next(g)
(2, 2, 1)

sum_of_three_squares

diofant.solvers.diophantine.sum_of_three_squares(n)[source]

Returns a 3-tuple \((a, b, c)\) such that \(a^2 + b^2 + c^2 = n\) and \(a, b, c \geq 0\).

Returns (None, None, None) if \(n = 4^a(8m + 7)\) for some \(a, m \in Z\). See [R482] for more details.

Parameters:

n : int

a non-negative integer.

References

[R482](1, 2) Representing a number as a sum of three squares, [online], Available: http://www.schorn.ch/howto.html

Examples

>>> from diofant.solvers.diophantine import sum_of_three_squares
>>> sum_of_three_squares(44542)
(207, 37, 18)

sum_of_four_squares

diofant.solvers.diophantine.sum_of_four_squares(n)[source]

Returns a 4-tuple \((a, b, c, d)\) such that \(a^2 + b^2 + c^2 + d^2 = n\).

Here \(a, b, c, d \geq 0\).

Parameters:

n : int

is a non-negative integer.

References

[R483]Representing a number as a sum of four squares, [online], Available: http://www.schorn.ch/howto.html

Examples

>>> from diofant.solvers.diophantine import sum_of_four_squares
>>> sum_of_four_squares(3456)
(8, 48, 32, 8)
>>> sum_of_four_squares(1294585930293)
(0, 1137796, 2161, 1234)

Internal Functions

These functions are intended for the internal use in Diophantine module.

merge_solution

diofant.solvers.diophantine.merge_solution(var, var_t, solution)[source]

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.

divisible

diofant.solvers.diophantine.divisible(a, b)[source]

Returns \(True\) if a is divisible by b and \(False\) otherwise.

extended_euclid

diofant.solvers.diophantine.extended_euclid(a, b)[source]

For given a, b returns a tuple containing integers \(x\), \(y\) and \(d\) such that \(ax + by = d\). Here \(d = gcd(a, b)\).

Parameters:

a, b : Integer

Returns:

tuple

\(x\), \(y\) and \(\gcd(a, b)\).

Examples

>>> from diofant.solvers.diophantine import extended_euclid
>>> extended_euclid(4, 6)
(-1, 1, 2)
>>> extended_euclid(3, 5)
(2, -1, 1)

PQa

diofant.solvers.diophantine.PQa(P_0, Q_0, D)[source]

Returns useful information needed to solve the Pell equation.

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 [R484] for more detailed information.

Parameters:

P_0, Q_0, D : Integer

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.

References

[R484](1, 2) Solving the generalized Pell equation x^2 - Dy^2 = N, John P. Robertson, July 31, 2004, Pages 4 - 8. http://www.jpr2718.org/pell.pdf

Examples

>>> from diofant.solvers.diophantine import PQa
>>> pqa = PQa(13, 4, 5)  # (13 + sqrt(5))/4
>>> next(pqa)  # (P_0, Q_0, a_0, A_0, B_0, G_0)
(13, 4, 3, 3, 1, -1)
>>> next(pqa)  # (P_1, Q_1, a_1, A_1, B_1, G_1)
(-1, 1, 1, 4, 1, 3)

equivalent

diofant.solvers.diophantine.equivalent(u, v, r, s, D, N)[source]

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 [R485]. No test is performed to test whether \((u, v)\) and \((r, s)\) are actually solutions to the equation. User should take care of this.

Parameters:u, v, r, s, D, N : Integer

References

[R485](1, 2) Solving the generalized Pell equation x**2 - D*y**2 = N, John P. Robertson, July 31, 2004, Page 12. http://www.jpr2718.org/pell.pdf

Examples

>>> from diofant.solvers.diophantine import equivalent
>>> equivalent(18, 5, -18, -5, 13, -1)
True
>>> equivalent(3, 1, -18, 393, 109, -4)
False

simplified

diofant.solvers.diophantine.simplified(x, y, z)[source]

Simplify the solution \((x, y, z)\).

parametrize_ternary_quadratic

diofant.solvers.diophantine.parametrize_ternary_quadratic(eq)[source]

Returns the parametrized general solution for the ternary quadratic equation eq which has the form \(ax^2 + by^2 + cz^2 + fxy + gyz + hxz = 0\).

References

[R486]The algorithmic resolution of Diophantine equations, Nigel P. Smart, London Mathematical Society Student Texts 41, Cambridge University Press, Cambridge, 1998.

Examples

>>> from diofant.abc import x, y, z
>>> from diofant.solvers.diophantine import parametrize_ternary_quadratic
>>> parametrize_ternary_quadratic(x**2 + y**2 - z**2)
(2*p*q, p**2 - q**2, p**2 + q**2)

Here \(p\) and \(q\) are two co-prime integers.

>>> parametrize_ternary_quadratic(3*x**2 + 2*y**2 - z**2 - 2*x*y + 5*y*z - 7*y*z)
(2*p**2 - 2*p*q - q**2, 2*p**2 + 2*p*q - q**2, 2*p**2 - 2*p*q + 3*q**2)
>>> parametrize_ternary_quadratic(124*x**2 - 30*y**2 - 7729*z**2)
(-1410*p**2 - 363263*q**2, 2700*p**2 + 30916*p*q - 695610*q**2, -60*p**2 + 5400*p*q + 15458*q**2)

diop_ternary_quadratic_normal

diofant.solvers.diophantine.diop_ternary_quadratic_normal(eq)[source]

Solves the quadratic ternary diophantine equation, \(ax^2 + by^2 + cz^2 = 0\).

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.

Examples

>>> from diofant.abc import x, y, z
>>> from diofant.solvers.diophantine import diop_ternary_quadratic_normal
>>> diop_ternary_quadratic_normal(x**2 + 3*y**2 - z**2)
(1, 0, 1)
>>> diop_ternary_quadratic_normal(4*x**2 + 5*y**2 - z**2)
(1, 0, 2)
>>> diop_ternary_quadratic_normal(34*x**2 - 3*y**2 - 301*z**2)
(4, 9, 1)

ldescent

diofant.solvers.diophantine.ldescent(A, B)[source]

Uses Lagrange’s method to find a non trivial solution to \(w^2 = Ax^2 + By^2\).

Here, \(A \neq 0\) and \(B \neq 0\) and \(A\) and \(B\) are square free. Output a tuple \((w_0, x_0, y_0)\) which is a solution to the above equation.

References

[R487]The algorithmic resolution of Diophantine equations, Nigel P. Smart, London Mathematical Society Student Texts 41, Cambridge University Press, Cambridge, 1998.
[R488]Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin, Mathematics of Computation, Volume 00, Number 0.

Examples

>>> from diofant.solvers.diophantine import ldescent
>>> ldescent(1, 1)  # w^2 = x^2 + y^2
(1, 1, 0)
>>> ldescent(4, -7)  # w^2 = 4x^2 - 7y^2
(2, -1, 0)

This means that \(x = -1, y = 0\) and \(w = 2\) is a solution to the equation \(w^2 = 4x^2 - 7y^2\)

>>> ldescent(5, -1)  # w^2 = 5x^2 - y^2
(2, 1, -1)

gaussian_reduce

diofant.solvers.diophantine.gaussian_reduce(w, a, b)[source]

Returns a reduced solution \((x, z)\) to the congruence \(X^2 - aZ^2 \equiv 0 \ (mod \ b)\) so that \(x^2 + |a|z^2\) is minimal.

Here w is a solution of the congruence \(x^2 \equiv a \ (mod \ b)\)

References

[R489]Gaussian lattice Reduction [online]. Available: http://home.ie.cuhk.edu.hk/~wkshum/wordpress/?p=404
[R490]Efficient Solution of Rational Conices, J. E. Cremona and D. Rusin, Mathematics of Computation, Volume 00, Number 0.

holzer

diofant.solvers.diophantine.holzer(x_0, y_0, z_0, a, b, c)[source]

Simplify the solution \((x_{0}, y_{0}, z_{0})\) of the equation \(ax^2 + by^2 = cz^2\) with \(a, b, c > 0\) and \(z_{0}^2 \geq \mid ab \mid\) to a new reduced solution \((x, y, z)\) such that \(z^2 \leq \mid ab \mid\).

prime_as_sum_of_two_squares

diofant.solvers.diophantine.prime_as_sum_of_two_squares(p)[source]

Represent a prime \(p\) which is congruent to 1 mod 4, as a sum of two squares.

References

[R491]Representing a number as a sum of four squares, [online], Available: http://www.schorn.ch/howto.html

Examples

>>> from diofant.solvers.diophantine import prime_as_sum_of_two_squares
>>> prime_as_sum_of_two_squares(5)
(2, 1)

pairwise_prime

diofant.solvers.diophantine.pairwise_prime(a, b, c)[source]

Transform \(ax^2 + by^2 + cz^2 = 0\) into an equivalent equation \(a'x^2 + b'y^2 + c'z^2 = 0\) where \(a', b', c'\) are pairwise relatively prime.

Returns a tuple containing \(a', b', c'\). \(\gcd(a, b, c)\) should equal \(1\) for this to work. The solutions for \(ax^2 + by^2 + cz^2 = 0\) can be recovered from the solutions of \(a'x^2 + b'y^2 + c'z^2 = 0\).

Examples

>>> from diofant.solvers.diophantine import pairwise_prime
>>> pairwise_prime(6, 15, 10)
(5, 2, 3)

make_prime

diofant.solvers.diophantine.make_prime(a, b, c)[source]

Transform the equation \(ax^2 + by^2 + cz^2 = 0\) to an equivalent equation \(a'x^2 + b'y^2 + c'z^2 = 0\) with \(\gcd(a', b') = 1\).

Returns a tuple \((a', b', c')\) which satisfies above conditions. Note that in the returned tuple \(\gcd(a', c')\) and \(\gcd(b', c')\) can take any value.

Examples

>>> from diofant.solvers.diophantine import make_prime
>>> make_prime(4, 2, 7)
(2, 1, 14)

reconstruct

diofant.solvers.diophantine.reconstruct(a, b, z)[source]

Reconstruct the \(z\) value of an equivalent solution of \(ax^2 + by^2 + cz^2\) from the \(z\) value of a solution of a transformed version of the above equation.

transformation_to_normal

diofant.solvers.diophantine.transformation_to_normal(eq)[source]

Returns the transformation Matrix from general ternary quadratic equation \(eq\) to normal form.

General form of the ternary quadratic equation is \(ax^2 + by^2 cz^2 + dxy + eyz + fxz\). This function returns a 3X3 transformation Matrix which transforms the former equation to the form \(ax^2 + by^2 + cz^2 = 0\). This is not used in solving ternary quadratics. Only implemented for the sake of completeness.