%I #30 Aug 20 2021 04:20:45
%S 3,19,163,571,1459,8803,9137,17497,41113,52489,78787,87211,135433,
%T 139483,144667,164617,174763,196579,274081,370009,370387,478243,
%U 760267,941489,944803,1041619,1220347,1236787,1319323,1465129,1663579,1994659
%N Prime factors of numbers in A006521 (numbers k that divide 2^k + 1).
%C A prime p is in this sequence iff all prime divisors of ord_p(2)/2 are in this sequence, where ord_p(2) is the order of 2 modulo p. - _Max Alekseyev_, Jul 30 2006
%H Joerg Arndt, <a href="/A057719/b057719.txt">Table of n, a(n) for n = 1..220</a> (terms up to 10^9, terms for n = 1..100 from T. D. Noe)
%H Alexander Kalmynin, <a href="https://arxiv.org/abs/1611.00417">On Novák numbers</a>, arXiv:1611.00417 [math.NT], 2016. See Chapter 4 p. 7 Novák primes.
%H C. Smyth, <a href="https://cs.uwaterloo.ca/journals/JIS/VOL13/Smyth/smyth2.html">The terms in Lucas Sequences divisible by their indices</a>, JIS 13 (2010) #10.2.4.
%e 2^171 + 1 == 0 (mod 171), 171 = 3^2*19, 2^13203+1 == 0 (mod 13203), 13203 = 3^4*163.
%t S = {2}; Reap[For[p = 3, p < 2 10^6, p = NextPrime[p], f = FactorInteger[ MultiplicativeOrder[2, p]]; If[f[[1, 1]] != 2 || f[[1, 2]] != 1, Continue[]]; f = f[[All, 1]]; If[Length[Intersection[S, f]] == Length[f], S = Union[S, {p}]; Print[p]; Sow[p]]]][[2, 1]] (* _Jean-François Alcover_, Nov 11 2018, from PARI *)
%o (PARI) { A057719() = local(S,f); S=Set([2]); forprime(p=3,10^7, f=factorint(znorder(Mod(2,p))); if(f[1,1]!=2||f[1,2]!=1,next); f=f[,1]; if(length(setintersect(S,Set(f)))==length(f), S=setunion(S,[p]); print1(p,", "))) }
%Y Cf. A006521, A066364.
%Y Cf. A136474, A136473.
%K nonn
%O 1,1
%A _Ignacio Larrosa Cañestro_, Oct 26 2000
%E Edited by _Max Alekseyev_, Jul 30 2006