# Maple coding to produce b-files for lists of # numbers with the same number of divisors. # Richard J. Mathar, http://www.strw.leidenuniv.nl/~mathar, 2007-03-24 # Return a list of of lists which are factorizations (product representations) # of n. Within each sublist, the factors are sorted. A minimum factor in # each element of sublists returned can be specified with 'mincomp'. # If mincomp=2, the number of sublists contained in the list returned is A001055(n). # Example: # n=8 and mincomp=2 return [[2,2,2],[4,8],[8]] listProdRep := proc(n,mincomp) local dvs,resul,f,i,j,rli,tmp ; resul := [] ; # list returned is empty if n < mincomp if n >= mincomp then if n = 1 then RETURN([1]) ; else # compute the divisors, and take each divisor # as a head element (minimum element) of one of the # sublists. Example: for n=8 use {1,2,4,8}, and consider # (for mincomp=2) sublists [2,...], [4,...] and [8]. dvs := numtheory[divisors](n) ; for i from 1 to nops(dvs) do # select the head element 'f' from the divisors f := op(i,dvs) ; # if this is already the maximum divisor n # itself, this head element is the last in # the sublist if f =n and f >= mincomp then resul := [op(resul),[f]] ; elif f >= mincomp then # if this is not the maximum element # n itself, produce all factorizations # of the remaining factor recursively. rli := listProdRep(n/f,f) ; # Prepend all the results produced # from the recursion with the head # element for the result. for j from 1 to nops(rli) do tmp := [f,op(op(j,rli))] ; resul := [op(resul),tmp] ; od ; fi ; od ; fi ; fi ; RETURN(resul) ; end: # Return a list of of lists which are ordered factorizations (product # representations) of n. The sublists, contain factors in all possible ways # of permutations. A minimum factor in each element of sublists returned can # be specified with 'mincomp'. # If mincomp=2, the number of sublists contained in the list returned is A07206(n). # Example: # n=8 and mincomp=2 return [[2,2,2],[4,8],[8,4],[8]] listProdRepSort := proc(n,mincomp) local usrt,resul,i,j,tmp ; # we first produce the list of factorizations with # elements sorted in the sublists. usrt := listProdRep(n,mincomp) ; resul := [] ; for i from 1 to nops(usrt) do # For each of the sorted sublists the result of the # case taking the order into account is produced by # taking all possible permutations (possibly only a single # one if all factors are the same). tmp := combinat[permute](op(i,usrt)) ; resul := [op(resul),op(tmp)] ; od ; RETURN(resul) ; end: # Produce a list of numbers with a prime power signature 'sig' # under the constraint that no number is larger than 'maxa' and # the minp'th prime is the smallest prime which is a factor of the # results. The signature 'sig' is a sorted list of numbers, each larger than zero, # with all prime power exponents of the resulting list of numbers. # Example: sig=[3,2,5] with minp=3 (indicating the 3rd prime, 5) tabulates # all numbers of the form p1^3*p2^2*p3^5 which are not larger than maxa, # where 5<=p1 maxa then break ; else # In a recursive construction, we drop the first # prime power (the first number) in the signature, # and look for all numbers that are residual factors # specified by the remaining exponents in the signature. if nops(sig) > 1 then # delete the first element for the sub signature sigsub := subsop(1=NULL,sig) ; # recursively call this function, where the # maxa is reduced by what already is considered # a factor (ff), and where the lowest prime # index of the residual factor is one higher # than the current one. subl := divsFixForsig(sigsub,maxa/ff,he+1) ; for i from 1 to nops(subl) do if op(i,subl)*ff <= maxa then resul := [op(resul),op(i,subl)*ff] ; fi ; od ; else # if the signature has only one element, this # is already the largest prime power. resul := [op(resul),ff] ; fi ; # we interpret the signature now again by assuming # the smallest prime number in the factorization # is the next prime after what was just considered. he := he+1 ; fi ; od ; RETURN(resul) ; end: # Produce a list of numbers with N divisors, all numbers # smaller or equal to 'maxa' # Note that in the special case of N being prime, one could # simply list all N-1st powers of consecutive primes; the algorithm # here is more generic and produces the same result in these cases. divsFix := proc(N,maxa) local preps,rec,sig,ele,resul,e ; # First generate lists of exponents which # factorize N. # Example: N=8 returns [2,2,2],[2,4],[4,2],[8] in 'preps' preps := listProdRepSort(N,2) ; # Generally, the number of divisors of some n is given # by numtheory[tau[(n) or by the product over (e_i-1) where # the e_i are the exponents in the prime number decomposition of # n. Above, we have decomposed N in all possible products # of this kind, so subtracting 1 leaves us with the # powers (prime exponents) of all possible numbers with # N divisors. # Example: N=8 with 'preps' as above is intepreted as signatures # [1,1,1],[1,3],[3,1] or [7] of prime power exponents # after subtracting 1 from each individual component. resul := [] ; # we loop over all these prime power signatures and accumulate # all results (which have no duplicates since each number has # only one unique signature of this type). for rec from 1 to nops(preps) do # from this signature of the product representation # of N we subtract 1 from each element. sig := op(rec,preps) ; for e from 1 to nops(sig) do ele := op(e,sig) ; sig := subsop(e=ele-1,sig) ; od ; # In the example, 'sig' would now be [1,1,1] # or [1,3] etc. resul := [op(resul),op(divsFixForsig(sig,maxa,1))] ; od ; RETURN(sort(resul)) ; end: # Auxiliary format to produce a b-file in the OEIS style # for a number N of divisors, written to the file 'file'. # The standard offset of 1 can be shifted by using a nonzero 'o'. bFile := proc(N,file,o) local amax,resul,i ; # We start under the mild assumption that we produce # all numbers below 50 with N divisors, and # increase this upper limit until we got # at least 400 lines in the b-file. amax := 50 ; resul := divsFix(N,amax) ; while nops(resul) < 400 do amax := 2*amax ; resul := divsFix(N,amax) ; od ; # 'resul' is now the sorted list, which is # printed with index. for i from 1 to nops(resul) do fprintf(file,"%d %d\n",o+i,op(i,resul)) ; od ; end: A000040 := proc(o) bFile(2,"b000040.txt",o) ; end: A001248 := proc(o) bFile(3,"b001248.txt",o) ; end: A030513 := proc(o) bFile(4,"b030513.txt",o) ; end: A030514 := proc(o) bFile(5,"b030514.txt",o) ; end: A030515 := proc(o) bFile(6,"b030515.txt",o) ; end: A030516 := proc(o) bFile(7,"b030516.txt",o) ; end: A030626 := proc(o) bFile(8,"b030626.txt",o) ; end: A030627 := proc(o) bFile(9,"b030627.txt",o) ; end: A030628 := proc(o) bFile(10,"b030628.txt",o) ; end: A030629 := proc(o) bFile(11,"b030629.txt",o) ; end: A030630 := proc(o) bFile(12,"b030630.txt",o) ; end: A030631 := proc(o) bFile(13,"b030631.txt",o) ; end: A030632 := proc(o) bFile(14,"b030632.txt",o) ; end: A030633 := proc(o) bFile(15,"b030633.txt",o) ; end: A030634 := proc(o) bFile(16,"b030634.txt",o) ; end: A030635 := proc(o) bFile(17,"b030635.txt",o) ; end: A030636 := proc(o) bFile(18,"b030636.txt",o) ; end: A030637 := proc(o) bFile(19,"b030637.txt",o) ; end: A030638 := proc(o) bFile(20,"b030638.txt",o) ; end: ########### call all of these... # no need for another prime number list... b-file exists # A000040(0) ; # this already exists with 5000 elements b-file # A001248(0) ; A030513(0) ; A030514(0) ; # offset 0 for no obvious reason A030515(-1) ; # offset 0 for no obvious reason A030516(-1) ; # offset 0 for no obvious reason A030626(-1) ; # offset 0 for no obvious reason A030627(-1) ; # there is some irregular leading number 1 in that sequence... fprintf("b030628.txt","0 1\n") ; A030628(0) ; A030629(0) ; # offset 0 for no obvious reason A030630(-1) ; # offset 0 for no obvious reason A030631(-1) ; # offset 0 for no obvious reason A030632(-1) ; # offset 0 for no obvious reason A030633(-1) ; # offset 0 for no obvious reason A030634(-1) ; # offset 0 for no obvious reason A030635(-1) ; # offset 0 for no obvious reason A030636(-1) ; # offset 0 for no obvious reason A030637(-1) ; # offset 0 for no obvious reason A030638(-1) ;