%I #53 Jun 11 2019 12:13:44
%S 1,2,3,6,10,21,30,42,78,110,210,330,390,546,903,930,1218,1806,1830,
%T 2310,2530,2730,4134,4290,6090,6162,6510,7590,9030,10230,12090,12246,
%U 12810,14910,15834,20130,20670,22110,23478,23790,28938,30030,30810,43134
%N Numbers n such that n divides the denominator of 2n-th Bernoulli number.
%C Numbers n such that the congruence k^(2n+1) == k (mod n) is true for 1<=k<=n. - _Michel Lagneau_, May 02 2012
%C In 2005, B. C. Kellner proved E. W. Weisstein's conjecture that denom(B_n) = n only if n = 1806. - _Jonathan Sondow_, Oct 14 2013.
%H Joerg Arndt, <a href="/A106741/b106741.txt">Table of n, a(n) for n = 1..594</a> (terms <= 10^8)
%H Bernd C. Kellner, <a href="http://bernoulli.org/~bk/denombneqn.pdf">The equation denom(B_n) = n has only one solution</a>, preprint 2005.
%H Victor Miller, <a href="https://listserv.nodak.edu/cgi-bin/wa.exe?A2=NMBRTHRY;9ede91fb.1205">Re: Q about a property of Bernoulli denominators</a>, NMBRTHRY list, May 5, 2012
%H Eric Weisstein's World of Mathematics, <a href="http://mathworld.wolfram.com/BernoulliNumber.html">Bernoulli Number</a>
%p for n from 1 to 10000 do:
%p m:=2*n+1: i:=1:
%p for k from 1 to n while(k &^ m mod n =k) do: i:=i+1: od:
%p if i=n then print(n) fi:
%p od: # _Michel Lagneau_, May 02 2012
%p A106741_list := proc(searchlimit) local isA106741, i;
%p isA106741 := proc(n)
%p numtheory[divisors](2*n);
%p map(i->i+1,%);
%p select(isprime,%);
%p mul(i,i=%) mod n = 0;
%p if % then n else NULL fi end:
%p seq(isA106741(i),i=1..searchlimit) end:
%p A106741_list(30000); # _Peter Luschny_, May 04 2012
%t okQ[n_] := AllTrue[Range[n], PowerMod[#, 2n+1, n] == Mod[#, n]&];
%t Reap[For[n = 1, n < 50000, n++, If[okQ[n], Sow[n]]]][[2, 1]] (* _Jean-François Alcover_, Jun 11 2019, after _Michel Lagneau_ *)
%o (PARI) is_A106741(n)=denominator(bernfrac(2*n))%n==0 \\ _Charles R Greathouse IV_, May 02 2012
%o (PARI){ for (n=1, 10^6, m = 2*n + 1; for (k=2, n, if ( Mod(k,n)^m != k, next(2) ); ); print1(n,", "); ); } /* _Joerg Arndt_, May 04 2012 */
%o (PARI) is_A106741(n)={ my(m=2*n+1); for(k=2, n, Mod(k, n)^m - k & return); 1} /* more than twice faster (in PARI 2.4.2) than with "if(...)" */ \\ _M. F. Hasler_, May 06 2012
%Y Cf. A002445, A014117.
%K nonn
%O 1,2
%A _Benoit Cloitre_, May 15 2005
%E Terms a(19)-a(29) from _Michel Lagneau_, May 02 2012
%E Terms >= 10230 by _Joerg Arndt_, May 04 2012