okQ[k_] := Module[{pp, ee}, {pp, ee} = FactorInteger[k] // Transpose; Max[ee] < 2 && Mod[Total[PowerMod[#, k, k]& /@ pp], k] == 0]; Reap[For[k = 6, k < 4*10^6, k++, If[CompositeQ[k], If[okQ[k], Print[k]; Sow[k] ] ] ] ][[2, 1]] (*
Jean-François Alcover, Feb 15 2018 *)