(* Warning: this program is just a verification of the existing data
and should not be used to extend the sequence beyond a(28) *)
BellMod[k_, m_] := Mod[Sum[Mod[StirlingS2[k, j], m], {j, 1, k}], m];
BellMod[k_, 1] := BellB[k];
period[nn_List] := Module[{lgmin=2, lgmax=5, nn1},
lg=If[Length[nn]<=lgmax, lgmin, lgmax];
nn1 = nn[[1;; lg]];
km=Length[nn]-lg;
Catch[Do[If[nn1==nn[[k;; k+lg-1]], Throw[k-1]];
If[k==km, Throw[0]], {k, 2, km}]]];
dd[n_] := SelectFirst[Table[{d, n/d},
{d, Divisors[n][[2;; -2]]}], GCD@@#==1&];
a[1]=1;
a[p_?PrimeQ] := a[p] = (p^p-1)/(p-1);
a[n_/; n>4 && dd[n]!={}] := With[{g = dd[n]}, LCM[a[g[[1]]], a[g[[2]]]]];
a[n_/; MemberQ[FactorInteger[n][[All, 2]], 1]] := a[n]=
With[{pp = Select[FactorInteger[n], #1[[2]] ==1 &][[All, 1]]},
a[n/Times@@pp]*Times@@a/@pp];
a[n_/; n>4 && GCD @@ FactorInteger[n][[All, 2]]>1] := a[n]=
With[{g=GCD @@ FactorInteger[n][[All, 2]]}, n^(1/g)*a[n^(1-1/g)]];
a[n_] := period[Table[BellMod[k, n], {k, 1, 28}]];
