function [ d, v, w, cmin, cmax ] = diophantine_2d_nonnegative ( a, b )

%*****************************************************************************80
%
%% diophantine_2d_nonnegative finds nonnegative diophantine solutions.
%
%  Discussion:
%
%     We are given a Diophantine equation 
%
%       a1 x1 + a2 x2 = b
%
%     for which the coefficients a are positive integers, and
%     the right hand side b is a nonnegative integer.
%
%     We are seeking all nonnegative integer solutions x.
%
%     This function analyzes the coefficient vector a, and returns
%     a value d, a vector v, and a vector 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 an arbitrary integer coefficient.
%
%  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): a particular solution for c = 0.
%
%    integer w(n): the general solution basis vector.
%
%    integer cmin, cmax: the range of coefficients c for which
%    the solution is nonnegative.
%
  d = NaN;
  v = [];
  w = [];
  cmin = 0;
  cmax = -1;

  if ( any ( a <= 0 ) )
    error ( 'diophantine_2d_nonnegative: a has a nonpositive entry.' )
  end

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

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

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

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

  d = i4vec_gcd ( a );
  if ( mod ( b, d ) ~= 0 )
    return
%   error ( 'diophantine_2d_nonnegative: b is not divisible by gcd ( a ).' );
  end 
%
%  Dimension the matrix.
%
  n = length ( a );
  if ( n ~= 2 )
    error ( 'diophantine_2d_nonnegative: a must have length 2.' );
  end

  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_2d_nonnegative - 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' );
    return
  end

  v = (A(1,2:np1))' * f;

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

  cmin = - Inf;
  cmax =   Inf;

  for i = 1 : n
    if ( w(i) < 0 )
      cmax = min ( cmax, floor ( - v(i) / w(i) ) );
    elseif ( 0 < w(i) )
      cmin = max ( cmin, ceil ( - v(i) / w(i) ) );
    end
  end

  return
end

