Number of Irreducible Polynomial Factors of a Polynomial in $\mathbb{F}_p[X]$ This question raises from the fact that over the integers, the prime omega $\Omega(N)$ function tells us the total number of prime factors of a number $N$. This is because the number $N$ factors as $N=\prod_{i=1}^{r}p_i^{e_i}$. Then $\Omega(N)= \sum_{i=1}^r e_i$.
Right now, my research conduces me to estimate the (average) number of irreducible factors of a monic polynomial over $\mathbb{F}_p$.
This is, take a polynomial $f(X) \in \mathbb{F}_p[X]$. Factor it over $\mathbb{F}_p$ as $f(X) = \prod_{i=1}^r f_i(X)^{e_i}$, where $f_i(X)$ is irreducible over $\mathbb{F}_p$. And finally obtain $\Omega(f):=\sum_{i=1}^re_i$.
My question is how I can estimate the number of irreducible factors of any polynomial $f(X)$ over $\mathbb{F}_p$.
 A: A possible approach is to analyze the number $N_{k,d}$ (throughout the answer, $p$ is considered fixed, and its primality is not used, so it can be a power of a prime as well) of monic polynomials $f$ over $\mathbb{F}_p$ of degree $d$ with $\Omega(f)=k$. We assume $\deg 1=\Omega(1)=0$ for the constant polynomial $1$, so that $N_{k,d}$ is defined for $k,d\geqslant 0$. Then, $N_{1,d}$ is the number of irreducible monic polynomials of degree $d$, and $\sum_{k=0}^d N_{k,d}=p^d$.
This is a road to various averages expressible. Say, $\Omega(f)$ averaged over monic $f$ with $\deg f=d$ is $$\overline\Omega(d):=p^{-d}\sum_{k=1}^d kN_{k,d}.$$ We denote by $M$ the set of all monic polynomials (over $\mathbb{F}_p$), and by $I$ the subset of $M$ consisting of all irreducible polynomials. Under our agreements above, we have $1\in M$ and $1\notin I$.
The analysis is done using generating functions: $$G_\Omega(s,t)=\sum_{k,d\geqslant 0}N_{k,d}s^k t^d=\sum_{f\in M}s^{\Omega(f)}t^{\deg f},\\G_I(t)=\sum_{d\geqslant 1}N_{1,d}t^d=\sum_{f\in I}t^{\deg f}=\frac{\partial}{\partial s}G_\Omega(s,t)\Bigg|_{s=0},\\G_{\overline\Omega}(t)=\sum_{d\geqslant 1}\overline\Omega(d)t^d=\frac{\partial}{\partial s}G_\Omega\left(s,\frac{t}{p}\right)\Bigg|_{s=1}.$$
The crucial fact (the unlabelled multiset construction of Flajolet-Sedgewick, applied) is $$G_\Omega(s,t)=\exp\left\{\sum_{n\geqslant 1}\frac{s^n}{n}G_I(t^n)\right\}$$ (TODO: write an elaborated appendix if needed). It also allows to find $G_I$, since $\sum_{k=0}^d N_{k,d}=p^d$ implies $G_\Omega(1,t)=\sum_{d\geqslant 0}(pt)^d=(1-pt)^{-1}$. So, with $\ell(z):=-\log(1-z)$, we obtain $$\sum_{n\geqslant 1}\frac1n G_I(t^n)=\ell(pt)\implies G_I(t)=\sum_{n\geqslant 1}\frac{\mu(n)}{n}\ell(pt^n)$$ by a variant of the Möbius inversion. Since we don't really need $G_\Omega$, let's compute directly $$G_{\overline\Omega}(pt)=\exp\big(\ell(pt)\big)\sum_{n\geqslant 1}G_I(t^n)=(1-pt)^{-1}\sum_{n\geqslant 1}\ell(pt^n)\sum_{m\,\mid\,n}\frac{\mu(m)}{m};$$ since $\sum_{m\mid n}\big(\mu(m)/m\big)=\varphi(n)/n$ using Euler's totient function, we find finally $$\overline\Omega(d)=\sum_{n=1}^d\frac{1}{np^n}\sum_{m\,\mid\,n}\varphi(m)p^{n/m}=\sum_{m=1}^d\frac{\varphi(m)}{m}\sum_{n=1}^{\lfloor d/m\rfloor}\frac{p^{(1-m)n}}{n}.$$
For large $p$, this is $\overline\Omega(d)=H_d+\dfrac{1}{2p}+\mathcal{O}(p^{-2})$ where $H_d=\displaystyle\sum_{n=1}^d\frac1n=\log d+\gamma+\ldots$

With $\omega$ in place of $\Omega$, things would go harder. We would use the powerset construction (instead of the multiset construction), but with another generating function in place of $G_I$: $$G_\omega(s,t)=\exp\left\{\sum_{n\geqslant 1}(-1)^{n-1}\frac{s^n}{n}\sum_{k\geqslant 1}G_I(t^{nk})\right\},$$ with somewhat laborous number-theoretic computations. We find $G_{\overline\omega}(pt)=A(t)B(t)$, where
\begin{align*}
A(t)&=\exp\sum_{n\geqslant 1}a_n\ell(pt^n),&B(t)&=\sum_{n\geqslant 1}b_n\ell(pt^n),\\a_n&=\sum_{\substack{a,b,c\geqslant 1\\abc=n}}\frac{(-1)^{a-1}}{a}\frac{\mu(c)}{c},&b_n&=\sum_{\substack{a,b,c\geqslant 1\\abc=n}}(-1)^{a-1}\frac{\mu(c)}{c}.
\end{align*}
Now both $a_n$ and $b_n$ can be found using Dirichlet series. For $a_n$, the result is simple: $$\sum_{n\geqslant 1}\frac{a_n}{n^s}=\sum_{a,b,c\geqslant 1}\frac{(-1)^{a-1}}{a^{1+s}}\frac{1}{b^s}\frac{\mu(c)}{c^{1+s}}=(1-2^{-s})\zeta(s)\implies a_n=\begin{cases}1,&n\text{ is odd}\\0,&n\text{ is even}\end{cases};\\\sum_{n\geqslant 1}\frac{b_n}{n^s}=(1-2^{1-s})\frac{\zeta^2(s)}{\zeta(1+s)}=\frac{(1-2^{1-s})(1-2^{-1-s})}{(1-2^{-s})^2}\prod_{p\neq 2}\frac{1-p^{-1-s}}{(1-p^{-s})^2},$$ which gives, for odd $n=p_1^{r_1}\cdots p_k^{r_k}$ and $r>0$, $$b_n=\prod_{j=1}^k\left[1+\left(1-\frac{1}{p_j}\right)r_j\right],\quad b_{2^r n}=-\frac{r}{2}b_n.$$
Thus, $G_{\overline\omega}(pt)=\left(\sum_{n\geqslant 1}b_n\ell(pt^n)\right)\prod_{n\geqslant 1}(1-pt^{2n-1})^{-1}$. This is where I got to at the moment.
