function [ d, v, W ] = diophantine_basis ( a, b )

%*****************************************************************************80
%
%% diophantine_basis() finds a basis for the solution space.
%
%  Discussion:
%
%     Given the Diophantine equation 
%
%       a1 x1 + a2 x2 + ... + an xn = b
%
%     we are seeking all possible integer solutions x.
%
%     This function analyzes the coefficient vector a, and returns
%     a value d, an n vector v, and an n x n-1 matrix W so that:
%     * The value b must be divisible by d;
%     * v is a solution;
%     * General solutions have the form x = v + W * c
%       where c is a set of n-1 arbitrary integer coefficients.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    11 March 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer a(n): the coefficients of the Diophantine equation.
%
%    integer b: the right hand side.
%
%  Output:
%
%    integer d: a value that must be a factor of the right hand side.
%
%    integer v(n): a particular solution for c = 0.
%
%    integer W(n,n-1): the basis vectors for the general solution.
%
  if ( any ( a == 0 ) )
    error ( 'diophantine_basis: a has a zero entry.' )
  end
%
%  Dimension the matrix.
%
  n = length ( a );
  np1 = n + 1;
%
%  Form the matrix: 
%
  A(1:n,1) = a(1:n);
  A(1:n,2:np1) = eye(n,n);
%
%  Perform repeated row reduction.
%
  while ( 1 < nnz ( A(1:n,1) ) )
%
%  Locate row containing nonzero entry of smallest magnitude.
%
    t = Inf;
    p = 0;
    for i = 1 : n
      s = abs ( A(i,1) );
      if ( s ~= 0 )
        if ( s < t )
          p = i;
        end
      end
    end
%
%  Swap row 1 and row p.
%
    T = A(1,1:np1);
    A(1,1:np1) = A(p,1:np1);
    A(p,1:np1) = T;
%
%  fix(r) rounds towards zero.
%
    for i = 2 : n
      s = fix ( A(i,1) / A(1,1) );
      A(i,1:np1) = A(i,1:np1) - s * A(1,1:np1);
    end

  end

  d = A(1,1);
  f = b / d;
  if ( f * d ~= b )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'diophantine_basis - Warning!\n' );
    fprintf ( 1, '  b = %d is not a multiple of d = %d\n', b, d );
    fprintf ( 1, '  There are no integer solutions to this equation.\n' );
  end
  v = (A(1,2:np1))' * f;
  W = (A(2:n,2:np1))';

  return
end

