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)= limlimits_{zto 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?

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