// Magma program to find, for n = 0..nMax, all numbers j such that tau(j^n) = j. [1]; // row for n=0 nMax:=23; for n in [1..nMax] do rowNTerms:=[1]; // every row (including row 0) includes "1" as a term // find all remaining terms of row n exps:=[1]; // exponents in primeSig of j; start w/primeSig of a single prime partProdTau:=[n+1,1]; // partial product array // (element 1 is the full product) of tau(j^n) tau:=partProdTau[1]; p:=2; while n mod p eq 0 do p:=NextPrime(p); end while; primesNotDivN:=[p]; // array: primes not dividing n partProdMinJ:=[p,1]; // partial prod array (elmnt 1 is full product) of min j while (tau ge partProdMinJ[1]) or (exps[1] gt 1) do if tau ge partProdMinJ[1] then // not out of bounds F:=Factorization(tau); if #F eq #exps then expsInTau:=[]; for j in [1..#F] do expsInTau[j]:=F[j][2]; end for; expsInTau:=Reverse(Sort(expsInTau)); if expsInTau eq exps then // term found rowNTerms[#rowNTerms+1]:=tau; end if; end if; jAdv:=1; // advance the 1st exponent by 1 else // out of bounds, but not done if exps[#exps] lt exps[1] - 1 then // need to find the lowest j such that exps[j] < exps[1] - 1 jAdv:=2; while exps[jAdv] ge exps[1] - 1 do jAdv+:=1; end while; else // extend arrays by 1 element p:=NextPrime(primesNotDivN[#primesNotDivN]); while n mod p eq 0 do p:=NextPrime(p); end while; exps[#exps+1]:=0; partProdTau[#partProdTau+1]:=1; primesNotDivN[#primesNotDivN+1]:=p; partProdMinJ[#partProdMinJ+1]:=1; jAdv:=#exps; end if; end if; // advance the jAdv-th exponent by 1 exps[jAdv]+:=1; partProdTau[jAdv]:=partProdTau[jAdv+1]*(exps[jAdv]*n+1); partProdMinJ[jAdv]:=partProdMinJ[jAdv+1]*primesNotDivN[jAdv]^exps[jAdv]; // reset the lower-indexed exponents to match the jAdv-th for j in [jAdv-1..1 by -1] do exps[j]:=exps[j+1]; partProdTau[j]:=partProdTau[j+1]*(exps[j]*n+1); partProdMinJ[j]:=partProdMinJ[j+1]*primesNotDivN[j]^exps[j]; end for; tau:=partProdTau[1]; end while; rowNTerms:=Sort(rowNTerms); rowNTerms; end for;