Linear diophantine equations

lindiop.rst

Linear diophantine equations

A general linear nonhomogeneous diophantine equation is of the form

a_0x_0 + a_1x_1 +\cdots+ 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_0, \ldots, x_n) is a solution with the right hand side d, then a solution for a multiple c = kd is given by (kx_0, \ldots, 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 \mathbf{Z}^{n+1} of all families (x_0, \dots, 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_0, \ldots, l_n) of the homogeneous equation where the l_i are linear combinations of n parameters t_0, \ldots, t_{n-1} with integer coefficients.

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

a_0x_0 + a_1x_1 + \cdots + a_nx_n = a_0x_0 + gy,

where g = \gcd(a_1, \ldots, a_n) and y is an integer. Then d = \gcd(a_0, g) and by means of the extended Euclidean algorithm we can find integers x and y such that

a_0x + gy = d.

Now, if a_1y_1 + \cdots + 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_0x + a_1y_1y + \cdots + a_ny_ny = d.

If (x_0, x_1, \ldots, x_n) is a solution of the homogeneous equation, then

a_0x_0 + gy = 0

for some integer y. This is equivalent to

ax_0 + by = 0,

where a = a_0/d and b = g/d are relatively prime. Hence x_0 = bt_0 for some integer t_0. Conversely, if this is true, then y = -at_0 satisfies the equation, and we obtain

a_0bt_0 - a(a_1y_1 + \cdots + a_ny_n)t_0 = 0,

so (bt_0, -ay_1t_0, \ldots, -ay_nt_0) is a solution of the homogeneous equation.

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

a_1x_1 + \ldots + a_nx_n = 0.

Hence, if (l_1, \ldots, l_n) is a general solution of this equation in terms of the parameters t_1, \ldots, t_{n-1},

(bt_0, l_1 - ay_1t_0, \ldots, l_n - ay_nt_0)

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)])