Find $\gcd(p^n-1,p^m+1)$.

We consider the case $p\geq 3$. For $p=2$ only 1 part is different and we make a remake when applicable. Let $d=\gcd(m,n)$ and $m=ud,n=vd$. Without loss of generality we also assume that $m\leq n$.

Theorem
We will show the following: If $u,v \not \equiv 0 \pmod 2$, then $$ \gcd(p^m+1,p^n+1) = p^d+1 $$ otherwise exactly one of $u,v$ is even and $$ \gcd(p^m+1,p^n+1) = 2 $$ (If $p=2$, then the latter $\gcd$ is $1$ instead.)

Euclidean algorithm
The proof follows a series of $\gcd$ reductions as follows:
(a) If $m=n$ or $m=0$: then $\gcd(p^m+1,p^n+1)= p^m+1$ and $2$ respectively. (Latter is $1$ if $p=2$.)
(b) If $2m \leq n$: Let $r=n-2m$, then $$ \begin{align*} \gcd(p^m+1,p^n+1) &= \gcd(p^m+1,(p^n+1) - p^{n-m}(p^m + 1) + p^{n-2m}(p^m+1)) \\ &=\gcd(p^m+1,p^{n-2m}+1) \\ &= \gcd(p^m+1,p^r+1)\\ &= \gcd(p^r+1,p^m+1) \end{align*} $$ Update $(m,n)=(r,m)$.
(c) Finally, since $m\leq n$, we have $m < n < 2m$. Let $r=2m-n > 0$, then $$ \begin{align*} \gcd(p^m+1,p^n+1) &= \gcd(p^m+1,p^{n+r}+p^r) \\ &= \gcd(p^m+1, p^{2m}+p^r) \\ &= \gcd(p^m+1, (p^{2m}+p^r) - p^m(p^m+1) + (p^m+1)) \\ &= \gcd(p^m+1,p^r+1)\\ &= \gcd(p^r+1,p^m+1) \end{align*} $$ Update $(m,n)=(r,m)$.
Note: We may check that $\gcd(u,v)=1$ is maintained for each iteration.

Observe that for (c), $r = 2m-n < 2m-m = m < n$ so both (b) and (c) strictly reduces $(m,n)$ to smaller values. Hence if we always have $m\neq n$ then $m=0$ at some point. This shows that the sequence terminates.

Proof
We start with $$ \gcd(p^m+1,p^n+1) = \gcd(p^{ud}+1,p^{vd}+1) $$ and apply the reductions iteratively. At each iteration that is not (a), since $r=n-2m$ or $r=2m-n$, we have $$ \gcd(m,r) = \gcd(m,\pm(n-2m)) = \gcd(m,n) = d $$ Hence after update of $(m,n)$ we still have the form $$ \gcd(p^{ud}+1,p^{vd}+1) $$
i.e. $\gcd(m,n)=d$ is maintained.

Moreover, the parity of $u,v$ is swapped but otherwise unchanged. We see this immediately for the new $v$ since it equals the old $u$. For the new $u=r/d$, we have $$ \begin{align*} r &= \pm(n-2m) =\pm(vd-2ud)\\ r/d &= \pm(v-2u) \equiv v \pmod 2 \end{align*} $$ So parity of new $u$ equals parity of old $v$.

Now we go back to the original case of $u,v \not\equiv 0\pmod 2$. Since both $u,v$ has odd parity, $u\neq 0$ and we cannot terminate at the case where $m=0$. This means we terminate at the case where $m=n$. Each update iteration maintains $\gcd(u,v)=1$ so this is only possible if $m=n=d$. Hence the result is $$\gcd(p^m+1,p^n+1)=\gcd(p^d+1,p^d+1)=p^d+1$$

On the other hand, $u\not \equiv v \pmod 2$, so we cannot have the case $m=n$. Hence we must terminate at $m=0$, which gives us $$\gcd(p^m+1,p^n+1)=\gcd(p^0+1,p^n+1)=2$$ (If $p=2$ then $p^n+1$ is odd so the $\gcd$ is 1.)

An alternate answer for $\gcd(p^m+1,p^n+1)$. The other case can be converted to this form via $$ \gcd(p^n-1,p^m+1) = \gcd(-(p^n-1)+p^n(p^m+1),p^m+1) = \gcd(p^{m+n}+1,p^m+1) $$ We remark that $p$ can be any integer $\geq 1$ and not necessarily prime.

Let $d=\gcd(m,n)$ and $m=ud,n=vd$. Denote $$ x = p^d, \;\;\; M=p^m+1 = x^u+1, \;\;\;N=p^n+1 = x^v+1 $$ In the following the key method we use is the evaluation of a telescoping series $$ \begin{align*} (1+x^k)\sum_{i=0}^{l}(-x^k)^i &= (1+x^k)-(x^k+x^{2k})+(x^{2k}+x^{3k}) + \cdots +(-x^k)^l(1+x^k)\\ &= 1 + x^k(-x^k)^l \end{align*} $$ i.e. only first and last term remains.

Case 1: exactly one of $u,v$ is even
We show that $D=\gcd(M,N)=1,2$ for even and odd $p$ respectively.
Let $A,B$ be $$ \begin{align*} A &= \sum_{i=0}^{v-1} (-x^u)^i, & B&= \sum_{i=0}^{u-1}(-x^v)^i \end{align*} $$ Then $$ \begin{align*} MA+NB &= \left((x^u+1)\sum_{i=0}^{v-1}(-x^u)^i\right) + \left((x^v+1)\sum_{i=0}^{u-1}(-x^v)^i \right) \\ &= (1+x^u(-x^u)^{v-1}) + (1+x^v(-x^v)^{u-1})\\ &= 2+(-1)^{v-1}x^{uv}+(-1)^{u-1}x^{uv}\\ &= 2 \end{align*} $$ since $u,v$ have different parity. Therefore $D=\gcd(M,N)$ must divide 2. If $p$ is odd then $M,N$ are both even, giving $D=2$. Otherwise if $p$ is even, $M,N$ are both odd so we must have $D=1$.

Case 2: $u,v$ both odd
We show that $D=\gcd(M,N) = p^d+1$ for any $p$.
Since $\gcd(u,v)=1$, we can find $r,s>0$ such that $$ 1+ur = vs $$ As both $u,v$ are odd, $$ 1+r \equiv s \pmod 2 $$ i.e. $r$ and $s$ have different parity. Now let $A,B$ be $$ \begin{align*} A &= \sum_{i=0}^{r-1}(-x^u)^i &, \;\;\;\;B &= \sum_{i=0}^{s-1}(-x^v)^i \end{align*} $$ Then $$ \begin{align*} xMA+NB &= x\left((x^u+1)\sum_{i=0}^{r-1}(-x^u)^i\right)+\left( (x^v+1)\sum_{i=0}^{s-1}(-x^v)^i\right)\\ &= x(1+x^u(-x^u)^{r-1}) + (1 + x^v(-x^v)^{s-1})\\ &= x + (-1)^{r-1}x^{ur+1} + 1 + (-1)^{s-1} x^{vs}\\ &= (x+1) + (-1)^{r-1}x^{vs} + (-1)^{s-1}x^{vs}\\ &= x+1 \\ &= p^d+1 \end{align*} $$ since $r,s$ have different parity. Once again $D=\gcd(M,N)$ must divide $p^d+1$. However since $m=ud,n=vd$ and $u,v$ are both odd, we can show that (*) $p^d+1 | p^{ud}+1$ and $p^d+1 | p^{vd}+1$. Therefore $D = p^d+1$. $$\tag*{$\Box$}$$

(*) Since $u-1$ is even, we have $$ \begin{align*} (1+x)\sum_{i=0}^{u-1} (-x)^i &= 1 + x(-x)^{u-1} = 1 + x^u \end{align*} $$ Same applies for $v$.