Writing a program that checks whether a linear equation has positive integer solutions

I'm trying to write an algorithm to determine whether a linear equation, specifically in the form of ax + by = c, has positive integer solutions for given a,b,c. It needs to be efficient since the numbers a, b, and c can be in range 0<=a,b,c<=10^16. How would I approach this problem?

Since it's a diophantine equation in question, I tried to check if the GCD of a and b divides c, but this way I can't differentiate between positive, negative, or zero solutions.

Algorithm to determine non-negative-values solution existance for linear diophantine equation

I found a solution here, but I didn't quite understand it. Maybe someone can simplify it for me? Since this one is pretty general and I'm only interested in equations with 2 variables.

Existence of integer solutions

Bezout's identity tells you indeed that with gcd(a,b) the greatest common divisor of a and b:

i) gcd(a,b) is the smallest positive integer that can be written as ax + by, and
ii) every integer of the form ax + by is a multiple of gcd(a,b).

So there you go. If c is dividable by gcd(a,b), you have solutions.

Finding positive solutions

From any pair of solutions, we can get all the other ones. Thus we can see if they can be positive. Still from the same identity, we get:

When one pair of Bézout coefficients (x0, y0) has been computed (e.g., using extended Euclidean algorithm), all pairs can be represented in the form

And now we're done. All you have to do is :

1. Use the extended Euclidean algorithm, that will give you
   • gcd(a,b) and
   • a pair of (x0,y0) such that a * x0 + b * y0 = gcd(a,b)
2. Check if gcd(a,b) divides c.
   • If not, no solutions exist.
   • If it does, multiply x0 and y0 by c / gcd(a,b) to get solutions for your equation.
3. If x0 and y0 both have the same sign, you're done.
   • If they are both positive, you have positive solutions,
   • If they are both negative, you don't.
4. If x0 and y0 have different signs, choose the smallest (in absolute value) k to change the sign of the negative integer.
   • That is, if x0 is negative, then take k = floor(d * x0 / b) (rounding to -infinity)
     and if is y0 is negative, take k = ceil(-d * y0 / a)
   • Compute (x1,y1) = (x0 - k * b / d , y0 + k * a / d)
     • If x1 and y1 are both positive, you just found two positive integer solutions.
     • If flipping the sign of one number flipped the other one, you can not find positive solutions.

Note that it is related to the question you linked, but the number of variables is different. This is solved because you only have two variables.