How can we compute the multiplicative partition function

There is a formula that depends only on the shape of the prime factorization of $n$, i.e. the exponents of the primes in the factorization. This formula is obtained with the Polya Enumeration Theorem (PET) applied to a certain repertoire. It includes the factorization into a product of one term, namely $n$.

Let $$n= p_1^{v_1} p_2^{v_2} \times \cdots \times p_q^{v_q}$$ be the prime factorization of $n$ and recall the function $\Omega(n)$ which is given by $$\Omega(n) = \sum_{k=1}^q v_k,$$ i.e. it counts the prime factors according to their multiplicities.

The repertoire that goes into PET is given by $$Q = -1 + \prod_{k=1}^q \sum_{v=0}^{v_k} u_{p_k}^v = -1 + \prod_{k=1}^q \frac{u_{p_k}^{v_k+1}-1}{u_{p_k}-1}$$ which is an ordinary generating function in the set of variables $\{u_{p_k}\}.$

With these settings the value $q(n)$ of the multiplicative partition function is given by $$q(n) = [u_{p_1}^{v_1} u_{p_2}^{v_2}\cdots u_{p_q}^{v_q}] \sum_{m=1}^{\Omega(n)} Z(S_m)(Q),$$ where $Z(S_m)$ is the cycle index of the symmetric group acting on $m$ elements, which is computed from $$Z(S_0) = 1 \quad \text{and}\quad Z(S_m) = \frac{1}{m} \sum_{l=1}^m a_l Z(S_{m-l}).$$

This produces the following sequence starting at $n=2$ and extending to $n=64$ $$1, 1, 2, 1, 2, 1, 3, 2, 2, 1, 4, 1, 2, 2, 5, 1, 4, 1, 4,\\ 2, 2, 1, 7, 2, 2, 3, 4, 1, 5, 1, 7, 2, 2, 2, 9, 1, 2, 2, 7,\\ 1, 5, 1, 4, 4, 2, 1, 12, 2, 4, 2, 4, 1, 7, 2, 7, 2, 2, 1, 11, 1, 2, 4, 11,\ldots$$ which points to OEIS A001055.

The complexity of the repertoire that goes into the cycle index is $\tau(n).$ Remark as of Tue Jan 7 00:40:27 CET 2014. This algorithm is definitely not as good as what was posted at the Math Overflow link. Please consult my second answer for a fast algorithm.

Here is the Maple code that was used to compute the above values:

with(numtheory);
with(group):
with(combinat):

pet_cycleind_symm :=
proc(n)
        option remember;

        if n=0 then return 1; fi;

        expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;

pet_varinto_cind :=
proc(poly, ind)
           local subs1, subs2, polyvars, indvars, v, pot, res;

           res := ind;

           polyvars := indets(poly);
           indvars := indets(ind);

           for v in indvars do
               pot := op(1, v);

               subs1 :=
               [seq(polyvars[k]=polyvars[k]^pot,
               k=1..nops(polyvars))];

               subs2 := [v=subs(subs1, poly)];

               res := subs(subs2, res);
           od;

           res;
end;

v :=
proc(n)
        option remember;
        local fact, rep, gf, fcount, p, res;

        fact := op(2, ifactors(n));

        rep := mul(add(u[fact[k][1]]^p, p=0..fact[k][2]),
                   k=1..nops(fact));
        rep := rep-1;

        res := 0;
        for fcount to bigomega(n) do
            gf := pet_varinto_cind(rep, pet_cycleind_symm(fcount));
            gf := expand(gf);

            for p to nops(fact) do
                gf := coeff(gf, u[fact[p][1]], fact[p][2]);
            od;

            res := res+gf;
        od;

        res;
end;