Proving that if $ad-bc = \pm 1$, then $\gcd(x,y) = \gcd(ax +by, cx + dy)$

Let $t = \gcd(x,y)$, then it is not hard to see that $t \mid ax +by$ and $t \mid cx+dy$ and hence $t \mid \gcd(ax +by, cx + dy)$.

Why does this work easily? Because we clearly see that $ax+by$ and $cx+dy$ can be obtained by taking a linear combination with coefficients in $\mathbb{Z}$.

If we could obtain $x$ and $y$ as a combination of $x' = ax+by$ and $y' =cx+dy$, then we would get the divisibility relation in the other direction in the same way.

Is this possible?

That is can we find integers $a', b', c', d'$ such that $$a' x' + b'y' = x$$ $$c' x' + d'y' = y$$

This is indeed guaranteed by the condition imposed, and can be checked by an explicit calculation as in another answer, so I will not repeat it, yet instead point out that if you know about matrices and determinants it is direct:

Namely setting $A = \begin{pmatrix}a & b \\ c & d \end{pmatrix}$, we have $$A\begin{pmatrix}x \\ y \end{pmatrix} = \begin{pmatrix}x' \\ y' \end{pmatrix}$$ and we seek a matrix $A'$ such that $$A'\begin{pmatrix}x' \\ y' \end{pmatrix} = \begin{pmatrix}x \\ y \end{pmatrix}$$ and this independent of the specific value of $x,y$. That means $A'A$ must be the identity. In other words we need $A$ to be invertible as an integer matrix, that is its inverse is also integral.

This is equivalent to its determinant being invertible as an integer, that is it is $\pm 1$, which is exactly the imposed condition.

This way of looking at it also shows how to generalized to more than two numbers.

$\forall m,n; \gcd(x,y)| (mx + ny)$

$\gcd(x,y)|(ax+by)$ and $\gcd(x,y)|(cx+dy)$

$\gcd(x,y)|\gcd(ax+by, cx+dy)$

$\gcd(ax+by, cx+dy)|(d(ax+by) - b(cx+dy))\\ \gcd(ax+by, cx+dy)|(ad - bc)x\\ \gcd(ax+by, cx+dy)|x$

similarly we can show that

$\gcd(ax+by, cx+dy)|(c(ax+by) - a(cx+dy))$ and $\gcd(ax+by, cx+dy)|y$

$\gcd(ax+by, cx+dy) | \gcd(x,y)$

if $\gcd(ax+by, cx+dy) | \gcd(x,y)$ and $\gcd(x,y)|\gcd(ax+by, cx+dy)$ then $\gcd(ax+by, cx+dy) = \gcd(x,y)$

It is the special case $\rm\, \Delta := {\rm det}\, A= \pm1\ $ in the following

Theorem $\ $ If $\rm\,(x,y)\overset{A}\mapsto (X,Y)\,$ is linear then $\: \rm\gcd(x,y)\mid \gcd(X,Y)\mid \Delta \gcd(x,y)$

Proof $\ $ Inverting the linear map $\rm\,A\,$ by Cramer's Rule (multiplying by the adjugate) yields

$$\rm \begin{eqnarray}\rm a\ x\, +\, b\ y &=&\rm X\\ \\ \rm c\ x\, +\, d\ y &\ =\ &\rm Y\end{eqnarray} \quad\Rightarrow\quad \begin{array} \rm\Delta\ x\ \ \ =\ \ \ \rm d\ X\, -\, b\ Y \\\\ \rm\Delta\ y\ =\ \rm -c\ X\, +\, a\ Y \end{array}\ ,\quad\ \Delta\ =\ ad-bc\qquad $$

Hence, by RHS system, $\rm\ n\ |\ X,Y\ \Rightarrow\ n\ |\ \Delta\:x,\:\Delta\:y\ \Rightarrow\ n\ |\ gcd(\Delta\:x,\Delta\:y)\ =\ \Delta\ gcd(x,y)\:.$
In particular $\rm\ n = \gcd(X,Y) \mid \Delta\, \gcd(x,y).\ $

Further, by LHS system, $\rm\,n\mid x,y\ \Rightarrow\ n\mid X,Y\ \Rightarrow\ n\mid\gcd(X,Y).$
In particular $\rm\ n = gcd(x,y)\mid \gcd(X,Y).\ \ \ $ QED

Remark $\ $ See here for many special-case applications.