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

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.