Is it true that $$\gcd\left(5^{2^n} + 1, 13^{2^n} + 1\right) = 2$$ for all $n \in \mathbf{Z}_{\geq 0}$?
I'm continually stumped with this and verifying it numerically is quite expensive very quickly (around $n = 20$).
I have previously blogged about the origins of the problem and included a proof of an easier version, but this continues to stump me.
Letting $F_n = 5^{2^n} + 1$ and $T_n = 13^{2^n} + 1$ it's clear both satisfy the recurrence: $$f(n) \left(f(n) - 2\right) = f(n + 1) - 2$$ but as of yet, I've been unable to utilize this. However, it may be useful in reducing the amount of computation for attacking it numerically. Concretely, starting from a known output from the Euclidean algorithm $$a_n F_n + b_n T_n = 2$$ the recurrence may be useful in determining $a_{n + 1}, b_{n + 1}$ in less expensive ways. (Of course we need to define these in such a way that they are unique, e.g. $0 \leq b_n < F_n$.)
[UPDATE]: I have verified @Zander's answer that both $F_{2206}$ and $T_{2206}$ are divisible by the value $p = 3 \cdot 2^{2208} + 1$ with the short snippet of Python below. This clearly shows the gcd is not $2$. I was also able to show that $p$ is prime (though this is not necessary) by using Proth's theorem and $a = 11$ as a witness to primality. Finally, I remain unclear on how such an $n$ (and with it a $p$) could be found.
n = 2206
modulus = 3 * (2 ** (n + 2)) + 1
f_residue = 5**(2**0) + 1
t_residue = 13**(2**0) + 1
for _ in xrange(n):
f_updated = f_residue * (f_residue - 2) + 2
f_residue = f_updated % modulus
t_updated = t_residue * (t_residue - 2) + 2
t_residue = t_updated % modulus
print '(5^(2^2206) + 1) MODULO (2^(2208) + 1):'
print f_residue
print '(13^(2^2206) + 1) MODULO (2^(2208) + 1):'
print t_residue
No, it's not true. $p=3\cdot 2^{2208}+1$ divides the GCD when $n=2206$.
If $p=3\cdot 2^{n+d}+1$ with $d>0$ then $x^{2^n}+1\equiv 0\pmod{p}$ has $2^n$ distinct roots in $[1,p-1]$. For $d=1$ this means $1/6$ of all residues are roots, and if you imagine their appearance in this class "looks random" in some sense, then for about $1/36$ of such primes both $5$ and $13$ will be roots.
Similarly when $d>1$ and also for primes of the form $k\cdot 2^{n+d}+1$ we can optimistically expect a fixed fraction of primes in the pattern to give counterexamples.
There's no guarantee because there may be some hidden structure, but I thought it likely that I could find a small-ish counterexample. Unfortunately I couldn't find one with $d=1$.