// Magma program to compute terms of OEIS A137825 nMax:=22; // # of terms to generate Prod:=1; maxMultiplicity:=0; SigmaPE:=[]; // SigmaPE[j] will store the sum of the divisors of // prime(j)^(e(j) + 1) where e(j) is the current exponent // of prime(j) in the prime factorization of Prod; P:=PrimesInInterval(2,NthPrime(nMax)); for j in [1..nMax] do SigmaPE[j]:=P[j]+1; end for; I:=[1]; // I[j] is the 1st prime index with multiplicity < j a:=[]; for n in [1..nMax] do currMultOpt:=0; jOpt:=I[currMultOpt+1]; SigmaMin:=SigmaPE[jOpt]; jtest:=jOpt; for currMult in [1..maxMultiplicity] do jtestPrev:=jtest; jtest:=I[currMult+1]; if jtest ne jtestPrev then if SigmaPE[jtest] le SigmaMin then // increasing currMult means decreasing jOpt, // so update even in case of a tie SigmaMin:=SigmaPE[jtest]; jOpt:=jtest; currMultOpt:=currMult; end if; end if; end for; if jOpt eq 1 then // one more 2 was included in Prod maxMultiplicity+:=1; I[maxMultiplicity+1]:=1; end if; I[currMultOpt+1]+:=1; SigmaPE[jOpt]:=SigmaPE[jOpt]*P[jOpt]+1; Prod:=Prod*P[jOpt]; a[n]:=Prod; end for; a;