function A = diophantine_basis2 ( a, b )

%*****************************************************************************80
%
%% diophantine_basis2() finds a basis for the solution space.
%
%  Discussion:
%
%     Given the Diophantine equations
%
%       a11 x1 + a12 x2 + ... + a1n xn = b1
%       ...
%       am1 x1 + am2 x2 + ... + amn xn = bm
%
%     we are seeking all possible integer solutions x.
%
%     This function analyzes the coefficient matrix a, and returns
%     a value d, an n vector v, and an n x n-1 matrix B so that:
%     * The value b must be divisible by d;
%     * v is a solution;
%     * General solutions have the form x = v + B * 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(m,n): the coefficients of the Diophantine equation.
%
%    integer b(m): the right hand side.
%
%  Output:
%
%    integer d: a value that must be a factor of the right hand side.
%
%    integer v(n): the solution for c = 0.
%
%    integer B(n,n-1): the basis vectors.
%

%
%  Dimension the matrix.
%  N = number of equations.
%  M = number of variables.
%
  [ n, m ] = size ( a );
%
%  Form the matrix [a';I]: 
%
  A(1:m,1:n) = a';
  A(1:m,n+1:n+m) = eye(m,m);
%
%  Perform repeated row reduction.
%
  for j = 1 : n

    while ( 1 < nnz ( A(j:m,j) ) )

      disp ( A )
      fprintf ( 'nnz = %d\n', nnz ( A(1:m,j) ) );
      pause ( )
%
%  Locate row containing nonzero entry of smallest magnitude.
%
      t = Inf;
      p = -1;
      for i = j : m
        s = abs ( A(i,j) );
        if ( s ~= 0 )
          if ( s < t )
            t = s;
            p = i;
          end
        end
      end
      fprintf ( ' Choose A(%d,%d)=%d\n', p, j, A(p,j) );
%
%  Swap row j and row p.
%
      T          = A(j,1:n+m);
      A(j,1:n+m) = A(p,1:n+m);
      A(p,1:n+m) = T;
%
%  fix(r) rounds towards zero.
%
      for i = j + 1 : m
        s = fix ( A(i,j) / A(j,j) );
        A(i,1:n+m) = A(i,1:n+m) - s * A(j,1:n+m);
      end

    end

  end

  return
end

