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));

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?

  • $\begingroup$ How do you feel about switch case statements? $\endgroup$ Jul 20 '20 at 18:55
  • $\begingroup$ @AbdullahAliSivas yes! that is an option I've been doing till now, Hard coding all the expressions of the residues. $\endgroup$
    – Galilean
    Jul 21 '20 at 1:07
  • $\begingroup$ 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 $\endgroup$ Jul 21 '20 at 23:42
  • $\begingroup$ 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\}$? $\endgroup$ Jul 23 '20 at 18:52
  • $\begingroup$ @FedericoPoloni for any other point which is not a pole, the residue function Will just output the number computed by the function. $\endgroup$
    – Galilean
    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.


Your Answer

