Of course, the first thing to do is find the prime factorisation of the number, like glowcoder said. Say
n = p^a * q^b * r^c * ...
Then
- find the multiplicative partitions of
m = n / p^a
- for
0 <= k <= a, find the multiplicative partitions of p^k, which is equivalent to finding the additive partitions of k
- for each multiplicative partition of
m, find all distinct ways to distribute a-k factors p among the factors
- combine results of 2. and 3.
It is convenient to treat the multiplicative partitions as lists (or sets) of (divisor, multiplicity) pairs to avoid producing duplicates.
I've written the code in Haskell because it's the most convenient and concise of the languages I know for this sort of thing:
module MultiPart (multiplicativePartitions) where
import Data.List (sort)
import Math.NumberTheory.Primes (factorise)
import Control.Arrow (first)
multiplicativePartitions :: Integer -> [[Integer]]
multiplicativePartitions n
| n < 1 = []
| n == 1 = [[]]
| otherwise = map ((>>= uncurry (flip replicate)) . sort) . pfPartitions $ factorise n
additivePartitions :: Int -> [[(Int,Int)]]
additivePartitions 0 = [[]]
additivePartitions n
| n < 0 = []
| otherwise = aParts n n
where
aParts :: Int -> Int -> [[(Int,Int)]]
aParts 0 _ = [[]]
aParts 1 m = [[(1,m)]]
aParts k m = withK ++ aParts (k-1) m
where
withK = do
let q = m `quot` k
j <- [q,q-1 .. 1]
[(k,j):prt | let r = m - j*k, prt <- aParts (min (k-1) r) r]
countedPartitions :: Int -> Int -> [[(Int,Int)]]
countedPartitions 0 count = [[(0,count)]]
countedPartitions quant count = cbParts quant quant count
where
prep _ 0 = id
prep m j = ((m,j):)
cbParts :: Int -> Int -> Int -> [[(Int,Int)]]
cbParts q 0 c
| q == 0 = if c == 0 then [[]] else [[(0,c)]]
| otherwise = error "Oops"
cbParts q 1 c
| c < q = [] -- should never happen
| c == q = [[(1,c)]]
| otherwise = [[(1,q),(0,c-q)]]
cbParts q m c = do
let lo = max 0 $ q - c*(m-1)
hi = q `quot` m
j <- [lo .. hi]
let r = q - j*m
m' = min (m-1) r
map (prep m j) $ cbParts r m' (c-j)
primePowerPartitions :: Integer -> Int -> [[(Integer,Int)]]
primePowerPartitions p e = map (map (first (p^))) $ additivePartitions e
distOne :: Integer -> Int -> Integer -> Int -> [[(Integer,Int)]]
distOne _ 0 d k = [[(d,k)]]
distOne p e d k = do
cap <- countedPartitions e k
return $ [(p^i*d,m) | (i,m) <- cap]
distribute :: Integer -> Int -> [(Integer,Int)] -> [[(Integer,Int)]]
distribute _ 0 xs = [xs]
distribute p e [(d,k)] = distOne p e d k
distribute p e ((d,k):dks) = do
j <- [0 .. e]
dps <- distOne p j d k
ys <- distribute p (e-j) dks
return $ dps ++ ys
distribute _ _ [] = []
pfPartitions :: [(Integer,Int)] -> [[(Integer,Int)]]
pfPartitions [] = [[]]
pfPartitions [(p,e)] = primePowerPartitions p e
pfPartitions ((p,e):pps) = do
cop <- pfPartitions pps
k <- [0 .. e]
ppp <- primePowerPartitions p k
mix <- distribute p (e-k) cop
return (ppp ++ mix)
It's not particularly optimised, but it does the job.
Some times and results:
Prelude MultiPart> length $ multiplicativePartitions $ 10^10
59521
(0.03 secs, 53535264 bytes)
Prelude MultiPart> length $ multiplicativePartitions $ 10^11
151958
(0.11 secs, 125850200 bytes)
Prelude MultiPart> length $ multiplicativePartitions $ 10^12
379693
(0.26 secs, 296844616 bytes)
Prelude MultiPart> length $ multiplicativePartitions $ product [2 .. 10]
70520
(0.07 secs, 72786128 bytes)
Prelude MultiPart> length $ multiplicativePartitions $ product [2 .. 11]
425240
(0.36 secs, 460094808 bytes)
Prelude MultiPart> length $ multiplicativePartitions $ product [2 .. 12]
2787810
(2.06 secs, 2572962320 bytes)
The 10^k are of course particularly easy because there are only two primes involved (but squarefree numbers are still easier), the factorials get slow earlier. I think by careful organisation of the order and choice of better data structures than lists, there's quite a bit to be gained (probably one should sort the prime factors by exponent, but I don't know whether one should start with the highest exponents or the lowest).
