First of all, if p divides a and b, it should also divide c, else equation has no solution.
We can then divide everything by p.
We care from now only of equations where a and b coprimes,
that is without any common divisor but 1.
Let ax_{0} + by_{0} = c one solution.
Then if ax + by = c is another one, by substracting comes :
a(x  x_{0}) + b(y  y_{0}) = 0. As a and b coprimes, a divides y  y_{0},
that is y  y_{0} = ka or :
y = y_{0} + ka then x = x_{0}  kb
Consider now equation ax + by = 1, still with a and b coprime.
The Bézout theorem says that two numbers a and b are coprime if and only if
there are two numbers u and v (in Z) with au + bv = 1.
This proves that equation ax + by = c, GCD(a, b) = 1
has always solutions.
If (u, v) is solution of ax + by = 1, then
(cu, cv) is solution of ax + by = c.
Summary :
Equation ax + by = c has solutions if and only if p = GCD(a,b) divides c
Then let a = pa', b = pb' c = pc', if x_{0},y_{0} is some solution of a'x + b'y = 1, then all solutions of ax + by = c are : x = c'x_{0}  kb' y = c'y_{0} + ka' for all k ∈ Z 
The problem is now to find some specific solution.
In very simple cases, we can just guess or try, but equations like
97x + 29y = 4000 already requires a more efficient method.
As 97 = 29×3 + 10 and 4000 = 29×137+27, equation 97x + 29y = 4000
can be written :
10x + 29(y + 3x  137) = 27 that is substituting z = y + 3x  137 : 10x + 29z = 27,
which is a simpler equation.
We can repeat in writing : 10(x + 2z  2) + 9z = 7 that is 10t + 9z = 7 then :
t + 9(z + t) = 7, or t + 9u = 7 which has obvious solution t = 7, u = 0
Doing variable substitution backwards, we get : z = u  t = 7,
then x = t + 2  2z = 23 and finally
y = z + 137  3x = 61, hence a specific solution x_{0} = 23, y_{0}=61,
and all solutions x = 23  29k, y = 61 + 97k.
x ≡ca^{φ(b)1} mod (b), y = (c  ax)/b 
Remains to calculate φ(b) (and for those who don't believe,
to prove the EulerFermat theorem)
A simple method to calculate φ(b) from the prime decomposition of b :
Let b = ∏p_{i}^{ai}, then
φ(b) = b∏(1  1/p_{i}) = ∏(p_{i}^{ai}  p_{i}^{ai1}).
Specifically, when p prime, φ(p) = p  1.
And for the above equation : x ≡ 4000×97^{φ(29)1} mod (29),
that is 4000×97^{27}
Hopefully, no need of calculate huge numbers as all calculations are modulo 29,
The numbers will never be more than 29² = 841.
4000 ≡ 27 mod (29) and 97 ≡ 10 mod (29)
already simplfy calculations into :
4000×97^{27} ≡ 27×10^{27} mod (29), Pfffiuu !
Remains to effectively calculate 10^{27} mod (29).
We use a method known as "Indian exponentiation"
Algorithm for calculation of V = A^{B}
Initialize V = 1
while B ≥ 1  if B odd, multiply V by A and decrease B by one  else, square A and divide B by 2 
Or a variante in binary :
Initialize V = 1
For each bit in exposant B, beginning from the most significant bit :  square V  if the bit is 1, multiply V by A 
which gives here 27 = 11011_{(2)} and 10^{27} = (((1²×10)²×10)²)²×10)²×10
As every multiplication or squaring is modulo 29 :
10^{2} = 100 ≡ 13 mod (29)
10^{3} ≡ 13×10 = 130 ≡ 14 mod (29)
10^{6} ≡ 14^{2} =196 ≡ 22 mod (29)
10^{12} ≡ 22^{2} =484 ≡ 20 mod (29)
10^{13} ≡ 20×10 = 200 ≡ 26 mod (29)
10^{26} ≡ 26^{2} = 676 ≡ 9 mod (29)
10^{27} ≡ 9×10 = 90 ≡ 3 mod (29)
And finally x ≡ 27×3 = 81 ≡ 23 mod (29)
Let : r_{0} = a, r_{1} = b (a > b > 0)
P_{0} = 0, P_{1} = 1 and Q_{0} = 1, Q_{1} = 0 Then repeat from i = 2
p = r_{n1} is GCD of a and b (Euclidean algorithm)

With 97x + 29y = 4000 this results into :
97 = 29×3 + 10, r_{2} = 10, q_{2} = 3, P_{2} = 3, Q_{2} = 1
29 = 10×2 + 9, r_{3} = 9, q_{3} = 2, P_{3} = 7, Q_{3} = 2
10 = 9×1 + 1, r_{4} = 1, q_{4} = 1, P_{4} = 10, Q_{4} = 3
9 = 1×9 + 0, r_{5} = 0, that's all and x = 3, y = 10 is solution of 97x + 29y = 1
Multiplying by c = 4000, which is a multiple of p = r_{4} = 1, we get
x = 12000 + 29k, y = 40000  97k
If we let m = k + [12000/29] = k + 413, we get back
x = 23 + 29m, y = 61  97m as above.
Script to solve ax + by = c with this algorithm.
x = 165  105m  7k
y = 110 + 70m + 5k z = 3m + 1 k and m independant in Z 
Let's give general formulas for solving equation ax + by + cz = D
Let p = GCD(a,b), a' = a/p, b' = b/p
Let u_{0} and v_{0} some solution of a'u + b'v = c z_{0}, t_{0} some solution of cz + pt = D x_{0}, y_{0} some solution of a'x + b'y = t_{0} The general solution of ax + by + cz = D is :

Direct application to 15x + 21y + 35z = 200 really gives :
p = 3, a' = 5, b' = 7 then 5u + 7v = 35 and u_{0} = 0 v_{0} = 5 (for instance)
35z + 3t = 200 z_{0} = 1, t_{0} = 55
5x + 7y = 55 : x_{0} = 11, y_{0} = 0 and finally :
x = 11 + 7k
y = 0  5k  5m z = 1 + 3m 
This method is quite tedious as it requires to successively solve several
equations in two unknowns, that is to perform the Euclide's algorithm several times.
It therefore becomes impracticable with more than 3 unknowns.
We then could efficiently solve this equation in considering it as a "system" (hum) of one équation (!)
in n unknowns...
(a  1 0) (b  0 1)Generalizing this reduction process to our equation 15x + 21y + 35z = 200, we perform a "row Euclidean reduction" of the matrix :
(15  1 0 0) (21  0 1 0) (35  0 0 1) resulting into : (15  1 0 0) (6 1 1 0) = minus [21/15]=1 times the 1st line (5 2 0 1) = minus [35/15]=2 times the 1st line (0  7 0 3) = minus [15/5]=3 times the 3rd line (1  1 1 1) = minus [6/5]=1 times the 3rd line (5 2 0 1) and finally : (0  7 0 3) (1  1 1 1) (0 7 5 6) = minus [5/5]=5 times the 2nd lineThe only line with a non null value in the 1st column (= GCD = 1) gives directly the Bézout relation : 15×1 + 21×1 + 35×(1) = PGCD = 1
The lines beginning with 0 in the above matrix give solutions of
the homogenous equation 15x + 21y + 35z = 0
That is (x = 7, y = 0, z = 3) and (x = 7, y = 5, z = 6), which are linearily independant,
and we may as well replace them by two independant linear combinations, for instance
u = (7, 5, 0) and v = (0, 5, 3)
We can then freely add any linear combination of u and v : ku + mv.
This gives the general solution of our equation :
x = 200 + 7k + 0m
y = 200  5k  5m z = 200 + 0k + 3m 
This method may be used for any number of unknowns.
Thus consider to write the system as matrix form (A)(X) = (B).
Perform the previous row Euclidean reduction on matrix
(A^{t}I), obtained by concatenation of the transposed of (A) and a unit matrix I.
We then get a matrix (RT), in which T is the transforming matrix : (RT) = (T)(A^{t}I)
(T) is invertible, as we just made linear combinations of rows.
This also results into (T)(A^{t}) = (R), hence :
(A)(T^{t}) = (R^{t}), and finally, as T is invertible :
(A) = (R^{t})(T^{t})^{1}
The determinant of T being ±1 (as product of matrices with determinant ±1, unimodular matrices),
the inverse (T^{t})^{1} is a matrix in integers.
Equation (A)(X) = (B) is then (R^{t})(T^{t})^{1}(X) = (B).
Let the auxiliary unknown (Y) = (T^{t})^{1}(X), this becomes
(R^{t})(Y) = (B), which is "simpler", because of the row Euclidean reduction of (R) into
a matrix which is at worse triangular.
Solving (A)(X) = (B) is then equivallent to solve (R^{t})(Y) = (B)
Then (Y) = (T^{t})^{1}(X) gives the general solution of (A)(X) = (B) :
(X) = (T^{t})(Y)
After this a little abstract proof of the method, just some example :
To solve the system
5x + 6y + 8z  2t = 3
6x  10y + 7z  8t = 11
That is in matrix form :
(5, 6, 8, 2)×(x) = (3 ) (6, 10, 7, 8) (y) (11) (z) (t)Perform the row Euclidean reduction of (A^{t}I) :
(5 6  1 0 0 0) (6 10  0 1 0 0) (8 7  0 0 1 0) (2 8  0 0 0 1) The lowest absolute value is 2 (1 10  1 0 0 2) (0 34  0 1 0 3) (0 25  0 0 1 4) (2 8  0 0 0 1) We may also change the signs of a row (1 10  1 0 0 2) (0 34  0 1 0 3) (0 25  0 0 1 4) (0 28 2 0 0 5)We thus get the GCD(5,6,8,2) = 1. Let's go on with the row reduction from the next columns in A^{t},
(1 10  1 0 0 2) (0 9  0 1 1 1) (0 25  0 0 1 4) (0 3 2 0 1 1) New min = 3 (1 1  5 0 3 1) We may also reduce this row to get a simpler R (diagonal) (0 0  6 1 2 4) (0 1  16 0 9 4) New min = 1 (0 3  2 0 1 1) (1 0  11 0 6 3 ) (0 0  6 1 2 4 ) (0 1  16 0 9 4 ) Or by exchanging rows : (0 0  50 0 28 13) (1 0  11 0 6 3 ) (0 1  16 0 9 4 ) (0 0  6 1 2 4 ) (0 0  50 0 28 13) (R) = (1 0) and (T) = (11 0 6 3 ) (0 1) (16 0 9 4 ) (0 0) (6 1 2 4 ) (0 0) (50 0 28 13)We have then to solve (R^{t})(Y) = (B) :
(1 0 0 0)×(y1) = (3 ) (0 1 0 0) (y2) (11) (y3) (y4)Which obviously gives y1 = 3 and y2 = 11, because R is "diagonal", if we can say that of a non square matrix.
(Y) = (3 ) (11) (k ) (m )and (X) = (T^{t})(Y) :
(x) = (11 16 6 50)×(3 ) = (209  6k  50m ) (y) (0 0 1 0 ) (11) (0 + k ) (z) (6 9 2 28) (k ) (117 + 2k + 28m) (t) (3 4 4 13) (m ) (53  4k  13m )
Of course, all these "filthy" calculations upon matrices are to be done
by a program,
as soon as the number of unknowns or of equations become a little large.
As in the above examples, the "simplest" writing for the solutions is not obvious.
Here, we may let m = 4 + m' :
x = 9  6k  50m'
y = k z = 5 + 2k + 28m' t = 1  4k  13m' 