Multiplication partitioning into k distinct elements

This problem is a straightforward application of the Polya Enumeration Theorem and very similar to the problem discussed at this MSE link.

Remark. What follows below does permit repeated factors, the formula for distinct factors is at the end.

Recall the recurrence by Lovasz for the cycle index $Z(S_n)$ of the multiset operator $\def\textsc#1{\dosc#1\csod} \def\dosc#1#2\csod{{\rm #1{\small #2}}}\textsc{MSET}_{=n}$ on $n$ slots, which is $$Z(S_n) = \frac{1}{n} \sum_{l=1}^n a_l Z(S_{n-l}) \quad\text{where}\quad Z(S_0) = 1.$$

We have for example, $$Z(S_3) = 1/6\,{a_{{1}}}^{3}+1/2\,a_{{2}}a_{{1}}+1/3\,a_{{3}}$$ and $$Z(S_4) = 1/24\,{a_{{1}}}^{4}+1/4\,a_{{2}}{a_{{1}}}^{2} +1/3\,a_{{3}}a_{{1}}+1/8\,{a_{{2}}}^{2}+1/4\,a_{{4}}.$$

Suppose the prime factorization of $n$ is given by $$n=\prod_p p^v.$$

Applying PET it now follows almost by inspection that the desired count is given by the term $$F(n, k) = \left[\prod_p X_p^v\right] Z(S_k)\left(\prod_p \frac{1}{1-X_p}\right)$$ where the square bracket denotes coefficient extraction of formal power series.

As an example of what these substituted cycle indices look like consider $$Z(S_3)\left(\frac{1}{1-X_1}\frac{1}{1-X_2}\frac{1}{1-X_3}\right) \\ = 1/6\,{\frac {1}{ \left( 1-X_{{1}} \right) ^{3} \left( 1-X_{{2}} \right) ^{3}\left( 1-X_{{3}} \right) ^{3}}} \\+1/2\,{\frac {1}{ \left( -{X_{{1}}}^{2}+1 \right) \left( -{X_{{2}}}^{2}+1 \right) \left( -{X_{{3}}}^{2}+1 \right) \left( 1-X_{{1}} \right) \left( 1-X_{{2}} \right) \left( 1-X_{{3}} \right) }} \\+1/3\,{\frac {1}{ \left( -{X_{{1}}}^{3}+1 \right) \left( -{X_{{2 }}}^{3}+1 \right) \left( -{X_{{3}}}^{3}+1 \right) }}.$$

It should be clear that coefficient extraction here is fast and efficient using the Newton binomial series which says that $$[Q^{\mu k}] \left(\frac{1}{1-Q^\mu}\right)^\nu = {k+\nu-1\choose \nu-1}.$$

We get for $F(n,1)$ the sequence $$1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,\ldots$$ which is as it ought to be.

For $F(n,2)$ we get the sequence $$1, 1, 1, 2, 1, 2, 1, 2, 2, 2, 1, 3, 1, 2, 2, 3, 1, 3, 1, 3,\ldots$$ which points us to OEIS A038548.

As another example consider $F(n,4)$ which yields $$1, 1, 1, 2, 1, 2, 1, 3, 2, 2, 1, 4, 1, 2, 2, 5, 1, 4, 1, 4,\ldots$$ which points us to OEIS A21320 where the above calculation is confirmed.

As a last example consider $F(n,5)$ which yields starting at $n=20$ $$4, 2, 2, 1, 7, 2, 2, 3, 4, 1, 5, 1, 7, 2, 2, 2, 9, 1, 2, 2,\ldots$$ which points us to OEIS A252665.

The following Maple code implements this cycle index computation.

with(combinat);
with(numtheory);

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;

F :=
proc(n, k)
    option remember;
    local v, q, sind, res;

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

    sind :=
    pet_varinto_cind(mul(1/(1-cat(`X`, q)), q=1..nops(v)),
                     pet_cycleind_symm(k));

    res := sind;
    for q to nops(v) do
        res := coeftayl(res, cat(`X`, q)=0, v[q][2]);
    od;

    res;
end;

Addendum. For distinct factors the formula is $$G(n, k) = \left[\prod_p X_p^v\right] Z(P_k)\left(\prod_p \frac{1}{1-X_p}\right)$$ where the square bracket denotes coefficient extraction of formal power series and $Z(P_n)$ is the cycle index of the set operator $\def\textsc#1{\dosc#1\csod} \def\dosc#1#2\csod{{\rm #1{\small #2}}}\textsc{SET}_{=n}$ which was also used in the linked-to computation from above.

Starting at $n=20$ we get for $G(n,3)$ the sequence $$2, 1, 1, 0, 4, 0, 1, 1, 2, 0, 4, 0, 2, 1, 1, 1, 4, 0, 1, 1,\ldots$$ which points us to OEIS A088434.

Starting at $n=120$ we get for $G(n,4)$ the sequence $$8, 0, 0, 0, 0, 0, 3, 0, 1, 0, 1, 0, 3, 0, 0, 1, 1, 0, 1, 0,\ldots$$

The Maple code for $G(n,k)$ is as follows.

with(combinat);
with(numtheory);

pet_cycleind_set :=
proc(n)
option remember;

    if n=0 then return 1; fi;

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

# duplicate code omitted

G :=
proc(n, k)
    option remember;
    local v, q, sind, res;

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

    sind :=
    pet_varinto_cind(mul(1/(1-cat(`X`, q)), q=1..nops(v)),
                     pet_cycleind_set(k));

    res := sind;
    for q to nops(v) do
        res := coeftayl(res, cat(`X`, q)=0, v[q][2]);
    od;

    res;
end;

Addendum. Observe that we can easily compute $$G(n,5)$$ where $$n=2^{50}\times 3^{20} \times 5^{10}$$ which has $11781$ divisors so that the brute force method would have to check $$1889560697901637761$$ five-tuples, which is not feasible. The answer is $$27879304592.$$