# Calculating residue of a rational function

I have a function $$f(z) = \frac{1}{(z-z_1)(z-z_2)(z-z_3)(z-z_4)}$$ All of $$\{z_1,z_2,z_3,z_4\}$$ are simple poles. The residues for this function are given as $$\text{Res}(f(z),z_i)= \lim\limits_{z\to zi} \frac{(z-z_i)}{(z-z_1)(z-z_2)(z-z_3)(z-z_4)}$$ For example to find $$\text{Res}(f(z),z_1)$$ one first cancels the $$(z-z_1)$$ numerator and denominators and then takes the limit $$z \to z_1$$. so the result is $$\text{Res}(f(z),z_1) = \frac{1}{(z_1-z_2)(z_1-z_3)(z_1-z_4)}$$ Similarly one can find $$\text{Res}(f(z),z_2)$$, $$\text{Res}(f(z),z_3)$$, $$\text{Res}(f(z),z_4)$$.

I want to implement this residue finding algorithm in a function. In Cpp I tried to implement this like this

double z1 = 1.0;
double z2 = 2.0;
double z3 = 3.0;
double z4 = 4.0;

auto res = [&](double z){
return [&](double zi){
return (z-zi)/((z-z1)*(z-z2)*(z-z3)*(z-z4));
}(z);
};


This returns -nan when I compute res(z1) as the function becomes of $$\frac{0}{0}$$ form. I wanted to define a function that will first get rid of the common factor in the numerator and the denominator and then puts the value $$z_1$$ in the function. For simple enough functions with simple poles, this should be enough to find the residue.

How to do this in Cpp?

• How do you feel about switch case statements? Jul 20 '20 at 18:55
• @AbdullahAliSivas yes! that is an option I've been doing till now, Hard coding all the expressions of the residues. Jul 21 '20 at 1:07
• That is pretty much the way to go. Doing symbolic computation is expensive and complicated, so if you have that solution already stick with it Jul 21 '20 at 23:42
• What is your input? The vector of the $z_i$s? And what is your output? If that function res is the proper signature, how is it supposed to behave for $zi \not \in \{z_1,z_2,z_3,z_4\}$? Jul 23 '20 at 18:52
• @FedericoPoloni for any other point which is not a pole, the residue function Will just output the number computed by the function. Jul 23 '20 at 20:08

Instead of hard coding all cases with a switch clause, you can parametrize the function by its poles:

double residue(size_t i, const std::vector<double> &poles) {
double res = 1.0;
for (size_t j=0; j < poles.size(); j++) {
if (j != i) {
res *= 1 / (poles[i] - poles[j]);
}
}
return res;
}


As a side note, I wonder whether de l'Hospital's rule might be helpful in case of less simple functions.