A partial fraction decomposition of an inverse of a generic polynomial of an arbitrary order. 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?
 A: 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:
In[1168]:= Clear[e]; Clear[a]; Clear[b]; Clear[j]; J =.; z =.; d = 4; \
MM = 2;
subst := j[0] :> 0;
Table[(Plus @@ 
     Table[Sum[((-1)^(e[i] - 1 - J - j[1])/(a[i] + b[i] z)^(1 + 
             j[1]) b[
            i]^(Plus @@ (e[#] & /@ Delete[Range[1, d], i]) - J) a[
            i]^(J + j[d] + 1 - e[i]) Binomial[J, 
           e[i] - 1 - 
            j[d]] Times @@ (Binomial[
                e[#[[2]]] - 1 + j[#[[1]] + 1] - j[#[[1]]], 
                e[#[[2]]] - 
                 1] b[#[[2]]]^(j[#[[1]] + 1] - 
                   j[#[[1]]])/(a[#[[2]]] b[i] - 
                   a[i] b[#[[2]]])^(e[#[[2]]] + j[#[[1]] + 1] - 
                   j[#[[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}] /. subst)]], {i, 1, d}] - z^J/
    Product[(a[xi] + z b[xi])^e[xi] , {xi, 1, d}]), 
  Evaluate[Sequence @@ 
    Join[Table[{e[xi], 1, MM}, {xi, 1, d}], {{J, 0, 
       d - 1}}]]] // Together


Out[1170]= {{{{{0, 0, 0, 0}, {0, 0, 0, 0}}, {{0, 0, 0, 0}, {0, 0, 0, 
     0}}}, {{{0, 0, 0, 0}, {0, 0, 0, 0}}, {{0, 0, 0, 0}, {0, 0, 0, 
     0}}}}, {{{{0, 0, 0, 0}, {0, 0, 0, 0}}, {{0, 0, 0, 0}, {0, 0, 0, 
     0}}}, {{{0, 0, 0, 0}, {0, 0, 0, 0}}, {{0, 0, 0, 0}, {0, 0, 0, 
     0}}}}}

