Let $d \ge 2 $ be an integer. Then let $\vec{n}:=\left(n_\xi \right)_{\xi=1}^d $ be integers (to be termed exponents) each of which is bigger or equal to one and then let $\left( a_\xi \right)_{\xi=1}^d $ and $\left( b_\xi \right)_{\xi=1}^d $ be parameters. We consider a polynomial whose roots are equal to $\left(-a_\xi/b_\xi \right)_{\xi=1}^d $ with multiplicities $\vec{n} $. Now, the partial fraction decomposition of the inverse of that polynomial reads:
\begin{eqnarray} &&\frac{1}{\prod\limits_{\xi=1}^d \left(a_\xi+b_\xi z \right)^{n_\xi}} = \\ &&\sum\limits_{i=1}^d \sum\limits_{0 \le j_1 \le j_2 \le \cdots \le j_{d-1} \le n_i-1} \left[ \prod\limits_{\begin{array}{lll} \xi=1,\xi\neq i \\ \lambda:= \xi - 1_{\xi \ge i+1}\end{array}}^d \binom{n_\xi+j_\lambda-j_{\lambda-1}-1}{n_\xi-1} \frac{b_\xi^{j_\lambda-j_{\lambda-1}}}{\left(a_\xi b_i - a_i b_\xi \right)^{n_\xi+j_\lambda-j_{\lambda-1}}} \right] \cdot \frac{(-1)^{j_{d-1}} b_i^{\sum\limits_{\xi=1,\xi\neq i}^d n_\xi}}{\left(a_i + b_i z \right)^{n_i-j_{j_{d-1}}}} \end{eqnarray}
We derived this formula by starting with the exponents being all equal to unity, then using the known partial fraction decomposition of that quantity, and then differentiating the result $(n_\xi-1)$-times with respect to $a_\xi$ and then dividing the outcome by $(-1)^{n_\xi-1} (n_\xi-1)! $. In all this $\xi=1,\cdots,d$. The code below verifies the correctness of this formula. We have:
In[232]:= Clear[j]; Clear[a]; Clear[b]; d = RandomInteger[{2, 4}]; z =.;
For[xi = 1, xi <= d, xi++, e[xi] = RandomInteger[{1, 5}]];
subst := j[0] :> 0;
res = Plus @@
Table[Sum[((-1)^(j[d - 1])/(a[i] + b[i] z)^(e[i] - j[d - 1]) b[
i]^(Plus @@ (e[#] & /@
Delete[Range[1, d], i])) Times @@ (Binomial[
e[#[[2]]] - 1 + j[#[[1]]] - j[#[[1]] - 1],
e[#[[2]]] -
1] b[#[[2]]]^(j[#[[1]]] -
j[#[[1]] - 1])/(a[#[[2]]] b[i] -
a[i] b[#[[2]]])^(e[#[[2]]] + j[#[[1]]] -
j[#[[1]] - 1]) & /@
Transpose[{Range[1, d - 1], Delete[Range[1, d], i]}]) /.
subst), Evaluate[
Sequence @@ (Table[{j[xi], j[xi - 1], e[i] - 1}, {xi, 1,
d - 1}] /. subst)]], {i, 1, d}];
Together[res] - Product[1/(a[xi] + b[xi] z)^e[xi], {xi, 1, d}]
Out[236]= 0
Having said all this my question would be the following. How would you go about generalizing this result to obtain a partial fraction decomposition of a generic rational function $P(z)/Q(z)$ where $P(z)$ and $Q(z)$ are polynomials of arbitrary orders?
This is not going to be the full answer to this question but only a part of it. Let us assume that the degree of the polynomial $P(z)$ in the numerator is strictly smaller than that of the one $Q(z)$ in the denominator. We also assume that we know the factorization of the numerator, i.e. $Q(z) = \prod\limits_{\xi=1}^d (a_\xi + b_\xi z)^{n_\xi} $. Then without loss of generality we can just take $P(z)= z^j $ where $j=0,\cdots, d-1 $. In this case the following decomposition holds:
\begin{eqnarray} &&\frac{z^j}{\prod\limits_{\xi=1}^d \left(a_\xi+b_\xi z \right)^{n_\xi}} = \\ &&\sum\limits_{i=1}^d \sum\limits_{0 \le j_1 \le j_2 \le \cdots \le j_{d} \le n_i-1} % \left[ \prod\limits_{\begin{array}{lll} \xi=1,\xi\neq i \\ \lambda:= \xi - 1_{\xi \ge i+1}\end{array}}^d \binom{n_\xi+j_{\lambda+1}-j_{\lambda}-1}{n_\xi-1} \frac{b_\xi^{j_{\lambda+1}-j_{\lambda}}}{\left(a_\xi b_i - a_i b_\xi \right)^{n_\xi+j_{\lambda+1}-j_{\lambda}}} \right] \cdot % (-1)^{n_i-1-j-j_1} b_i^{-j+\sum\limits_{\xi=1,\xi\neq i}^d n_\xi} a_i^{j+1+j_d-n_i} \binom{j}{n_i-1-j_d} \frac{1}{\left(a_i + b_i z \right)^{1+j_1}} \end{eqnarray}
We derived this result by starting with all the exponents $(n_i)_{i=1}^d $ equal to unity then using the known partial fraction decomposition in that case and then differentiating the both sides by $(n_\xi-1)$ times by $a_\xi$ and dividing the result by $(-1)^{n_\xi-1} (n_\xi-1)! $. All this for $\xi=1,\cdots,d$.
As always the script below verifies this formula: