Before posting my last version, I had tried with **ListTools:-BinarySearch**, but that consumed about half of the time. This time, I implemented a binary search tailored to this algorithm. As Joe said, it's probably not worth it; but it's probably not a significant waste either.

A much more significant change is that I made it call **numtheory:-divisors **only once and **sort **only once. All subsequent sorted sublists of divisors are made my filtering the this one list. And the filtered lists are stored in a local remember table.

It seems that there is considerable interest in this algorithm. It certainly fascinates me. It amazes me that simply keeping the lists of factorings carefully sorted and only splitting the last element is sufficient to generate all the factorings. Joe: Is your method for generating all partitions of a multiset similar to this?

**Factorings:= proc(n::posint, m::posint:= 0)**

** local **

** L, K, Div, i, j, k, Nsq, cnt**

** ,Omega:= numtheory:-bigomega(n) # max # of factors**

** ,DIVS:= sort([(numtheory:-divisors(n) minus {1,n})[]]) # not changed during this proc's run.**

** # Get sorted list of divisors of divisors by filtering DIVS; also returns pointer to midpt of list.**

** ,divisors:= proc(n::posint)**

** option remember;**

** local D:= select(x-> irem(n,x)=0, DIVS)[1..-2]; # remove n itself from end.**

** D, iquo(nops(D)+1, 2) # Extra 1 is for square root.**

** end proc**

** # Find first index k into Div[1..Nsq] such that Div[k] >= x. **

** ,BinarySearch:= proc(x)**

** local lo, hi, m;**

** if Nsq=0 then return 1 elif Div[Nsq] < x then return Nsq+1 fi; # fall-thru cases. **

** lo:= 1; hi:= Nsq;**

** do**

** m:= iquo(hi+lo, 2);**

** if Div[m] = x then return m elif x < Div[m] then hi:= m-1 else lo:= m+1 fi;**

** if lo > hi then return lo fi **

** end do**

** end proc **

** ; **

** if m > Omega then return [] elif m > 0 then Omega:= m fi; **

** L[1]:= [[n]]; # Initialize loop w/ sole factoring of length 1.**

** # Avoid wasted effort of filtering divisor list on the first (i=2) pass.**

** divisors(n):= DIVS, iquo(nops(DIVS)+1, 2);**

** for i from 2 to Omega do # i is number of factors.**

** cnt:= 0; # cnt is index for list-building table.**

** # j is a single factoring for previous i, so j is a list.**

** for j in L[i-1] do**

** cnt:= cnt+1;**

** # Split last, and greatest, member of list, if not prime (prime case falls thru).**

** Div, Nsq:= divisors(j[-1]);**

** # j[-1] = Div[k]*Div[-k] for any k, since Div is sorted.**

** # Only use splits that maintain the list in order. **

** K[cnt]:= seq([j[1..-2][], Div[k], Div[-k]], k= `if`(i=2, 1, BinarySearch(j[-2]))..Nsq)**

** end do;**

** # Convert table K to list L[i]. **

** L[i]:= [seq](K[j], j= 1..cnt)**

** end do;**

** # Convert table L to list for final output.**

** `if`(m = 0, [seq](L[i][], i= 1..Omega), L[m])**

**end proc;**