Stable numerical partial fraction decomposition for large scale problem

Question

Given a rational polynomial $\frac{f(x)}{g(x)}$ with $f$ and $g$ having same degree and real coefficients. The degree of $f$ and $g$ is a few thousand. I would like to compute the partial fraction decomposition: $$ \frac{f(x)}{g(x)} = \sum_i \frac{r_i}{x - p_i}, $$ where $p_i$ and $r_i$ are the poles and residues of the rational polynomial.

Unfortunately, MATLABs residuez breaks down numerically in the degree of few hundred. Is there any way to perform the partial fraction decomposition numerically more stable? What is the proper way to assess the difficulty of the problem?

I am aware that the partial fraction decomposition is an ill-posed problem. Some more information on my problem: The poles lie on the unit circle and have multiplicity of 1. The poles can be accurately computed for a degree of few thousand. I only need the absolute value of the residues. None of the zeros and poles cancel each other. All coefficients lie between -1 and 1, and most of the coefficients are zero.

Answer 1

If these are simple poles, the residues are $$r_j = \lim_{x \to p_j} \frac{(x - p_j) f(x)}{g(x)} = \frac{f(p_j)}{g'(p_j)}$$ If the coefficients of $f$ and $g$ (and thus of $g'$) are large, accurate computation of $f(p_j)$ and $g'(p_j)$ will be difficult. I suggest you use a Computer Algebra System such as Maple or Mathematica that can use higher precision.