%I #18 May 24 2023 11:16:28
%S 7381,13357,18421,69121,70021,129961,192589,241501,271921,313501,
%T 342793,399421,423613,625681,780913,1265641,1362097,1501489,1566181,
%U 1673101,1691521,1728001,2228941,2381401,2472301,2642221,3156661,3171961,3383281,3557521,3730321,4033861,4233061,4831201,5387041,6720961
%N Composite numbers k such that b^k == b (mod sigma(k)) for every integer b.
%p N:= 10^7: # to get all terms <= N
%p Q:= select(isprime, [seq(q,q=5..N/9, 4)]):
%p filter:= proc(n) local s; uses numtheory;
%p s:= sigma(n);
%p issqrfree(s) and andmap(p -> (n-1 mod (p-1) = 0), factorset(s));
%p end proc:
%p select(filter, {seq(seq(q*m^2, m = 3..floor(sqrt(N/q)),2),q=[1,op(Q)])}); # _Robert Israel_, Sep 21 2016
%t M = 10^7; (* to get all terms <= M *)
%t Q = Select[Range[5, M/9, 4], PrimeQ];
%t filter[n_] := Module[{s = DivisorSigma[1, n]}, SquareFreeQ[s] && AllTrue[FactorInteger[s][[All, 1]], Mod[n-1, #-1] == 0&]];
%t Select[Union@Flatten[Table[Table[q*m^2, {m, 3, Floor[Sqrt[M/q]], 2}], {q, Prepend[Q, 1]}]], filter] (* _Jean-François Alcover_, May 24 2023, after _Robert Israel_ *)
%Y Cf. A002808, A000203.
%K nonn
%O 1,1
%A _Altug Alkan_, Oct 09 2016