Solution (source code)

= Solution
{c}

First recall why the two-variable <Bézout's identity> follows from the <Euclidean algorithm>. The initial values $a,b$ are <integer linear combinations> of $a,b$. At every subsequent step, $r_{j+1}=r_{j-1}-q_jr_j$ is again an <integer linear combination> of $a,b$. Therefore the last nonzero <remainder>, which is their <greatest common divisor>, has that form.

Put $e=\gcd(a,b)$. The preceding argument supplies <integers> $u,v$ with $e=au+bv$. Applying the same argument to $e,c$ gives <integers> $p,q$ such that
$$
\gcd(e,c)=ep+cq.
$$
An <integer> divides $a,b,c$ exactly when it divides $e,c$: if it divides $a,b$, it divides $e=au+bv$, while $e$ itself divides both $a,b$. Hence $\gcd(e,c)=\gcd(a,b,c)=d$. Substituting the expression for $e$ proves
$$
d=(au+bv)p+cq=a(up)+b(vp)+cq.
$$
Thus
$$
\boxed{x=up,\qquad y=vp,\qquad z=q}
$$
are <integers> satisfying $ax+by+cz=d$. This also explains how to construct the coefficients using two applications of the <Euclidean algorithm>.