function [ d, v, B, kmin, kmax ] = diophantine_nonnegative ( a, b )

%*****************************************************************************80
%
%% diophantine_nonnegative() finds nonnegative diophantine solutions.
%
%  Discussion:
%
%     We are given a Diophantine equation 
%
%       a1 x1 + a2 x2 + ... + an xn = b
%
%     for which the coefficients a and right hand side b are nonnegative,
%     and for which at least one coefficient is strictly positive.
%
%     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 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:
%
%    07 April 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): the solution for c = 0.
%
%    integer B(n,n-1): the basis vectors.
%
  if ( any ( a < 0 ) )
    error ( 'diophantine_nonnegative: a has a negative entry.' )
  end

  if ( ~ i4vec_is_integer ( a ) )
    error ( 'diophantine_nonnegative: some a is not an integer.' )
  end

  if ( sum ( a ) <= 0 )
    error ( 'diophantine_nonnegative: a does not have a positive entry.' )
  end

  if ( b < 0 )
    error ( 'diophantine_nonnegative: b is negative.' )
  end

  if ( ~ i4_is_integer ( b ) )
    error ( 'diophantine_nonnegative: b is not an integer.' )
  end

  d = i4vec_gcd ( a );
  if ( mod ( b, d ) ~= 0 )
    error ( 'diophantine_nonnegative: b is not divisible by gcd ( a ).' )
  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_solve - 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;

  B = (A(2:n,2:np1))';

  kmin = - Inf * ones ( n - 1, 1 );
  kmax =   Inf * ones ( n - 1, 1 );

  for j = 1 : n - 1
    for i = 1 : n
      if ( B(i,j) < 0 )
        kmax(j) = min ( kmax(j), - v(i) / B(i,j) );
      elseif ( 0 < B(i,j) )
        kmin(j) = max ( kmin(j), - v(i) / B(i,j) );
    end
  end

  return
end

