Linear diophantine equations

lindiop.rst

Linear diophantine equations

A general linear nonhomogeneous diophantine equation is of the form

a₀x₀ + a₁x₁ + ⋯ + a_nx_n = c

where the coefficients a_i and the constant c are integers. The components x_i of its solutions are also expected to be integers.

Solutions

The possible values of the left hand side of the equation form a subgroup of the additive group of integers. It consists of the multiples of its least positive element d, which is the greatest common divisor of the coefficients a_i.

Hence the equation is solvable if and only if c is divisible by d. Moreover, if (x₀, …, x_n) is a solution with the right hand side d, then a solution for a multiple c = kd is given by (kx₀, …, kx_n).

The difference of two solutions of the equation is a solution of the homogeneous diophantine equation whose right hand side is zero. Conversely, the sum of a particular solution of the equation and a solution of the homogeneous equation is a solution of the nonhomogeneous equation.

The solutions of the homogeneous equation form a subgroup of the group Z^n + 1 of all families (x₀, …, x_n) of integers. Assuming the coefficients a_i non-zero, as we may, the group of solutions is a free group of rank n. Hence its elements can be represented by n parameters taking integral values.

Implementation

For the complete solution of the nonhomogeneous linear diophantine equation it suffices to find

• the greatest common divisor d of the coefficients a_i,
• a particular solution of the equation with the right hand side d,
• a general solution (l₀, …, l_n) of the homogeneous equation where the l_i are linear combinations of n parameters t₀, …, t_n − 1 with integer coefficients.

These can be computed recursively as follows. For any integers x_i we have

a₀x₀ + a₁x₁ + ⋯ + a_nx_n = a₀x₀ + gy,

where g = gcd (a₁, …, a_n) and y is an integer. Then d = gcd (a₀, g) and by means of the extended Euclidean algorithm we can find integers x and y such that

a₀x + gy = d.

Now, if a₁y₁ + ⋯ + a_ny_n = g is a particular solution of the equation with n variables, we obtain a new particular solution for n + 1 variables

a₀x + a₁y₁y + ⋯ + a_ny_ny = d.

If (x₀, x₁, …, x_n) is a solution of the homogeneous equation, then

a₀x₀ + gy = 0

for some integer y. This is equivalent to

ax₀ + by = 0,

where a = a₀/d and b = g/d are relatively prime. Hence x₀ = bt₀ for some integer t₀. Conversely, if this is true, then y =  − at₀ satisfies the equation, and we obtain

a₀bt₀ − a(a₁y₁ + ⋯ + a_ny_n)t₀ = 0,

so (bt₀,  − ay₁t₀, …,  − ay_nt₀) is a solution of the homogeneous equation.

The difference of any two solutions of the homogeneous equation having the same first component bt₀ is a solution of the homogeneneous equation of n variables

a₁x₁ + … + a_nx_n = 0.

Hence, if (l₁, …, l_n) is a general solution of this equation in terms of the parameters t₁, …, t_n − 1,

(bt₀, l₁ − ay₁t₀, …, l_n − ay_nt₀)

is a general solution of the homogeneous equation of n + 1 variables.

This recursive procedure can be coded as follows.

def _linear_diop(coeffs, params):
  """Helper function for _diop_linear.

  Recursively solve a linear diophantine equation.

  Parameters
  ==========
  coeffs: Coefficients of a homogeneous linear equation.
          List of n + 1 non-zero integers.
  params: Parameters for the general homogeneous solution.
          List of n symbols.

  Returns
  =======
  d : int
      GCD of the coefficients.
  part : list of int
      Particular solution of the equation with right hand side ``d``.
  gen : list of Expr
      General solution of the homogeneous equation in terms of
      ``params``.
  """
  if len(coeffs) == 1:
      return abs(coeffs[0]), [sign(coeffs[0])], [0]

  g, part, gen = _linear_diop(coeffs[1:], params[1:])
  x, y, d = extended_euclid(coeffs[0], g)
  a, b = coeffs[0]//d, g//d
  t = params[0]

  return (d, [x] + [c*y for c in part],
          [b*t] + [-a*c*t + f for c, f in zip(part, gen)])

SymPy's _diop_linear may then be given the following form.

def _diop_linear(var, coeff, param):
  if not var:
      return ()

  coeffs = [coeff[v] for v in var]
  params = list(symbols(param.name +':' + str(len(var) - 1), integer=True))
  d, part, gen = _linear_diop(coeffs, params)

  if Integer(1) not in coeff:
      return tuple(gen)

  c = -coeff[Integer(1)]
  if c%d != 0:
      return tuple([None]*len(var))
  c //= d
  return tuple([c*x + f for x, f in zip(part, gen)])