$4x^2+1$ factors only into $4y+1$ primes How can one prove that numbers of the form $4x^2+1$ can only be divided by primes of the form $4y+1$ (e.g. there is no $x$ for which $7$ divides $4x^2+1$)? 
On a quick lookup, the statement is given as the first comment in https://oeis.org/A053755 and of course, manual testing concurs with this
 A: In general; note that your question is about the form $4x^2 + y^2$ with $\Delta = -16.$ If $q \equiv 3 \pmod 4$ is a prime that divides the form, then it divides both $x,y.$ However, in your case $y=1.$
Given a binary quadratic form
$$ f(x,y) = a x^2 + b xy+ c y^2   $$
with $a,b,c$ integers.
Given 
$$  \Delta = b^2 - 4 a c,  $$
where we require that $\Delta,$ if non-negative, is not a square (so also $\Delta \neq 0,1$).
Proposition: given an odd prime $q$ such that $q$ does not divide $\Delta$ and, in fact, $$ (\Delta| q) = -1, $$
 whenever $$ f(x,y) \equiv 0 \pmod q,   $$ then BOTH
$$ x \equiv 0 \pmod q, \; \; \; \; \; y \equiv 0 \pmod q.   $$
Proof: the integers taken $\pmod q$ form a field; every nonzero element has a multiplicative inverse. If either $a,c$ were divisible by $q,$ we would have $\Delta $ equivalent to $b^2$ mod $q,$ which would cause $(\Delta|q)$ to be $1$ rather than $-1.$ As a result, neither $a$ nor $c$ is divisible by $q.$
Next, $q$ is odd, so that $2$ and $4$ have multiplicative inverses mod $q,$ they are non zero in the field. Put togethaer, $$4a \neq 0 \pmod q$$
Now, complete the square:
$$ a x^2 + b xy+ c y^2 \equiv 0 \pmod q  $$
if and only if
$$ 4a(a x^2 + b xy+ c y^2) \equiv 0 \pmod q,  $$
$$ 4a^2 x^2 + 4abxy + 4ac y^2   \equiv 0 \pmod q,       $$
$$ 4a^2 x^2 + 4abxy + (b^2 y^2 - b^2 y^2) + 4ac y^2   \equiv 0 \pmod q,       $$
$$ 4a^2 x^2 + 4abxy + b^2 y^2 - (b^2 y^2 - 4ac y^2)   \equiv 0 \pmod q,       $$
$$ (4a^2 x^2 + 4abxy + b^2 y^2) - (b^2  - 4ac) y^2   \equiv 0 \pmod q,       $$
$$ (2ax + by)^2 - (b^2  - 4ac) y^2   \equiv 0 \pmod q,       $$
$$ (2ax + by)^2 - \Delta y^2   \equiv 0 \pmod q,       $$
$$ (2ax + by)^2 \equiv \Delta y^2   \pmod q.       $$
Now, ASSUME that $y \neq 0 \pmod q.$ Then $y$ has a multiplicative inverse which we are allowed to call $1/y,$ and we have
$$ \frac{(2ax + by)^2}{y^2} \equiv \Delta   \pmod q,       $$ 
$$ \left( \frac{2ax + by}{y} \right)^2 \equiv \Delta   \pmod q.       $$
However, the HYPOTHESIS that $ (\Delta| q) = -1 $ says that this is impossible, thus contradicting the assumption that $y \neq 0 \pmod q.$
So, in fact, $y \equiv 0 \pmod q.$ The original equation now reads
$$ a x^2  \equiv 0 \pmod q  $$ with the knowledge that $a \neq 0 \pmod q,$ so we also get
$$ x \equiv 0 \pmod q, \; \; \; \; \; y \equiv 0 \pmod q.   $$
$$ \bigcirc   \bigcirc  \bigcirc  \bigcirc  \bigcirc  \bigcirc  \bigcirc  \bigcirc $$
