How many ways to write a number $n$ as the product of natural numbers $\geq 2$? I am looking for a closed form (or efficient algorithm) for $f(n)$, the number of ways in which $n$ can be written as a product of natural numbers $\geq 2$. Note that $f(n)=\sum_{i=1}^{\Omega(n)}{g(n,i)}$, where $g(n,i)$ is the number of ways in which $n$ can be written as a product of $i$ natural numbers $\geq 2$ and $\Omega(n)$ is the total number of prime factors of $n$. If helpful, we can assume that we're given the prime factorization of $n$.
For example, if $n=40$, then we have the following:
$\Omega(40)=\Omega(2^35^1)=3+1=4$
$g(40,1)=1: (40)$
$g(40,2)=3: (20, 2), (10, 4), (8, 5)$
$g(40,3)=2: (10, 2, 2), (5, 4, 2)$
$g(40,4)=1: (5, 2, 2, 2)$
so $f(40)=1+3+2+1=7$
Partition (number theory) - Restricted part size or number of parts has potentially useful results. It's easy to get $g(n,2)$, as it's just $\frac{d(n)-2}{2}$, where $d(n)$ is the number of divisors of $n$.
 A: Among  several possible  approaches  we  can view  this  problem as  a
straightforward application of the  Polya Enumeration Theorem and very
similar     to    the     problem    discussed     at     this    MSE
link.
Recall  the recurrence by Lovasz  for the cycle  index $Z(S_n)$ of
the multiset operator $\mathfrak{M}_{=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.$$
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 of  multiplicative partitions into  $k$ factors is given  by the
term
$$F(n, k) = \left[\prod_p X_p^v\right]
Z(S_k)\left(-1 + \prod_p \frac{1}{1-X_p}\right)$$
where the square bracket denotes coefficient extraction 
of formal power series.
We are interested in
$$G(n) = \sum_{k=1}^{\Omega(n)} F(n, k).$$
Now recall the generating function of the cycle index of the symmetric
group which is
$$H(z) = \sum_{q\ge 0} Z(S_q) z^q =
\exp\left(a_1z + a_2\frac{z^2}{2} + a_3\frac{z^3}{3}+\cdots\right)
= \exp\left(\sum_{l\ge 1} a_l \frac{z^l}{l}\right).$$
This gives for $G(n)$
$$G(n) = \left[\prod_p X_p^v\right]
\exp\left(\sum_{l\ge 1} \frac{1}{l}
\left(-1 + \prod_p \frac{1}{1-X_p^l}\right)\right).$$
Recalling $v_p(n)$ which  is the exponent of the  maximum power of $p$
that divides $n$ we finally obtain
$$G(n) = \left[\prod_p X_p^v\right]
\exp\left(\sum_{l=1}^{\max_p v_p(n)}  \frac{1}{l}
\left(-1 + \prod_p \frac{1}{1-X_p^l}\right)\right).$$
This formula  has a certain  aesthetic value. Implemented in  Maple it
will produce the following sequence:
$$1, 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, \ldots$$
which  points  us  to OEIS  A001055  where
additional  information  awaits,   including  simple  recurrences  for
practical computation of this function.
The  Maple code  for the  above formula  is as  follows (here  we have
replaced the  prime number index on  the variable $X$  by a positional
index on the variable $Y$):

with(numtheory);

facts :=
proc(n)
    option remember;
    local f, gf, p, res;

    f := op(2, ifactors(n));
    gf := exp(add(1/l*
                  (-1+mul(1/(1-Y[p]^l), p=1..nops(f))),
                  l=1..max(map(p->p[2], f))));

    res := gf;
    for p to nops(f) do
        res := coeftayl(res, Y[p]=0, f[p][2]);
    od;

    res;
end;

