login
Composite numbers k such that b^k == b (mod sigma(k)) for every integer b.
0

%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