function x = diophantine_3d_positive ( a, b )

%*****************************************************************************80
%
%% diophantine_3d_positive finds positive diophantine solutions.
%
%  Discussion:
%
%     We are given a Diophantine equation 
%
%       a1 x1 + a2 x2 + a3 x3 = b
%
%     for which the coefficients a are positive integers, and
%     the right hand side b is a nonnegative integer.
%
%     We are seeking all strictly positive integer solutions x.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 May 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer a(3): the coefficients of the Diophantine equation.
%
%    integer b: the right hand side.
%
%  Output:
%
%    integer x(3,k): solutions to the equation.
%
  if ( length ( a ) ~= 3 )
    error ( 'diophantine_3d_nonnegative: a does not have length 3.' )
  end

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

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

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

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

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

  d = i4vec_gcd ( a );
  if ( mod ( b, d ) ~= 0 )
    error ( 'diophantine_3d_nonnegative: b is not divisible by gcd ( a ).' );
  end
%
%  Modify problem.
%
  b2 = b - sum ( a(1:3) );

  x = [];
%
%  Create reduced 2d equation.
%
  gr = gcd ( a(1), a(2) );
  ar = a(1:2) / gr;
  br = 1;

  as = [ 1; a(3) ];
  bs = b2;

  [ ds, vs, ws, cmins, cmaxs ] = diophantine_2d_nonnegative ( as, bs );

  for cs = cmins : cmaxs
    x3 = vs(2) + cs * ws(2);
    ar = a(1:2);
    br = b2 - a(3) * x3;
    if ( mod ( br, gr ) == 0 )
      [ dr, vr, wr, cminr, cmaxr ] = diophantine_2d_nonnegative ( ar, br );
      for cr = cminr : cmaxr
        x1 = vr(1) + cr * wr(1);
        x2 = vr(2) + cr * wr(2);
        x = [ x, [x1+1;x2+1;x3+1;] ];
      end
    end
  end

  return
end

