|
%I #16 Oct 26 2023 09:36:42
%S 2,3,4,5,6,12,20,30,72,165,288,693,1056,3024,9280,22500,42845,60480,
%T 240000,794580,1814400,7040040,26352000,98654400,321552000,1260230400,
%U 5311834416,17570520000,75087810000,325180275840,1526817600000
%N a(1) = 2; for n>1, a(n)=SENSigma(a(n-1)), where SENSigma(m) = (-1)^((Sum_i r_i)+Omega(m))*Sum_{d|m} (-1)^((Sum_j Max(r_j))+Omega(d))*d = Product_i (Sum_{1<=s_i<=r_i} p_i^s_i)+(-1)^(r_i+1) if m=Product_i p_i^r_i, d=Product_j p_j^r_j, p_j^max(r_j) is the largest power of p_j dividing m.
%C By "Max(r_j)" is meant the following: if d|m, d=p^e*q^f, m=p^x*q^y*r^z then Max(e)=x, Max(f)=y.
%p SENSigma := proc(n) local ifs,i,a,r,p ; ifs := ifactors(n)[2] ; a := 1 ; for i from 1 to nops(ifs) do r := op(2,op(i,ifs)) ; p := op(1,op(i,ifs)) ; a := a*(p*(1-p^r)/(1-p)-(-1)^r) ; od ; RETURN(a) ; end: A125141 := proc(nmax) local a ; a := [2] ; while nops(a)< nmax do a := [op(a),SENSigma(op(-1,a))] ; od ; RETURN(a) ; end: A125141(40) ; # _R. J. Mathar_, May 18 2007
%t SENSigma[n_] := Module[{Ifs, i, a, r, p }, Ifs = FactorInteger[n]; a = 1; For[i = 1, i <= Length[Ifs], i++, r = Ifs[[i, 2]]; p = Ifs[[i, 1]]; a = a(p(1 - p^r)/(1 - p) - (-1)^r)]; Return[a]];
%t A125141[nmax_] := Module[{a}, a = {2}; While[Length[a] < nmax, a = Append[a, SENSigma[a[[-1]]]]]; Return[a]];
%t A125141[40] (* _Jean-François Alcover_, Oct 26 2023, after _R. J. Mathar_ *)
%Y Cf. A126851, A126852, A125142.
%K nonn
%O 1,1
%A _Yasutoshi Kohmoto_, Jan 12 2007
%E Edited by _N. J. A. Sloane_ at the suggestion of _Andrew S. Plewe_, May 14 2007
%E More terms from _R. J. Mathar_, May 18 2007
|
