OFFSET
1,1
COMMENTS
For each prime p == 1 (mod 4), there are two values of x (mod p^5) that solve x^2 + 1 == 0 (mod p^5), and then x + k*p^5 is in the sequence for every k. Thus the asymptotic density of this sequence should be 1 - Product_p (1 - 2/p^5), where the product is over all primes p == 1 (mod 4). - Robert Israel, Sep 04 2018
LINKS
Robert Israel, Table of n, a(n) for n = 1..10000
EXAMPLE
1068 is in the sequence because 1068^2+1 = 1140625 = 5^6*73;
143044 is in the sequence because 143044^2+1 = 20461585937 = 13^5*55109;
390112 is in the sequence because 390112^2+1 = 152187372545 = 5*13*17^6*97.
MAPLE
N:= 10^5: # to get all terms <= N
P:= select(isprime, [seq(i, i=5..floor((N^2+1)^(1/5)), 4)]):
g:= proc(x, r, N) local t; t:= rhs(op(x)); seq(t+r*k, k=0..(N-t)/r) end proc:
R:= `union`(seq(map(g, {msolve(n^2+1, p^5)}, p^5, N), p=P)):
sort(convert(R, list)); # Robert Israel, Sep 04 2018
MATHEMATICA
Select[Range[2, 20000], Max[Transpose[FactorInteger[#^2+1]][[2]]]>4&]
PROG
(PARI) isok(n) = vecmax(factor(n^2+1)[, 2]) >= 5; \\ Michel Marcus, Sep 04 2018
CROSSREFS
KEYWORD
nonn
AUTHOR
Michel Lagneau, Nov 02 2012
STATUS
approved