?

No sé como se demuestra, pero lo estoy comprobando y de momento hasta n=150000 lo cumplen todos.

Bueno, tras tres horas voy a dejarlo, ha llegado hasta n=800000 y todos los n^2+1 cumplen que todos sus divisores primos impares son de la forma 4k+1. La probabilidad de que sea cierto al tenor de estos resultados es prácticamente total. He revisado el programa para ver si le puedo dar más velocidad, creo que calculaba muchas veces raíces cuadradas de forma innecesaria y no le puse optimizaciones al compilador.

Está bien la respuesta dada por Gregorio Morales pero poco explicada. El teorema de Euler dice que si n y p son coprimos entonces:

nφ(p)≡1(modp)nφ(p)≡1(modp)

Lo primero que debemos comprobar para poder aplicarlo es que n y p son coprimos. En la situación que plantea tenemos que p es un primo de la forma 4k+3, luego los únicos números que no serían coprimos con p tendrían p como factor primo.

Es decir, un número n no coprimo con p sería:

n=pj⋅qn=pj·q

con lo cual

n2+1=p2jq2+1≡0+1=1(modp)n2+1=p2jq2+1≡0+1=1(modp)

y esto es absurdo ya que por hipótesis se cumple

n2+1≡0(modp)n2+1≡0(modp)

luego es cierto que n y p serán coprimos con lo cual podemos asegurar que

nφ(p)≡1(modp)nφ(p)≡1(modp)

Aparte sabemos que

n4≡1(modp)n4≡1(modp)

pero eso por si solo no quiere decir que

4|φ(p)4|φ(p)

también podría ser que el indicador de Euler de p fuera múltiplo de 2 sin serlo de 4

Ahora bién, eso no puede ser ya que

n2≡−1(modp)n2≡−1(modp)

para que esta congruencia valga 1 debe elevarse el n2n2 a una potencia par 2m2m con lo cual

(n2)2m=n4m(n2)2m=n4m

con lo cual será

φ(p)=4mφ(p)=4m

al ser p primo

p−1=4mp−1=4m

p=4m+1p=4m+1