Consider the following integral in two complex variables $z_1$ and $z_2$:
$$\frac{1}{(2\pi i)^2}\oint_{{|z_1|=\epsilon}\atop{|z_2|=\epsilon}}dz_1 dz_2\frac{1}{z_1 z_2\left(1+a+\frac{z_1}{z_1+z_2}\right)\left(1+b+\frac{z_2}{z_1+z_2}\right)}$$
If I perform the contour integration in $z_1$ first, get the residue in $z_1$, and then perform the contour integration in $z_2$ to get the final residue, I get overall:
$$Res_{z_1=0,z_2=0}=\frac{1}{(1+a)(2+b)}$$
If I reverse the process and integrate the $z_2$ contour first, and then $z_1$, I end up with the final result:
$$Res_{z_1=0,z_2=0}=\frac{1}{(2+a)(1+b)}$$
These two overall residues are not the same, which implies that something must have gone wrong in the above calculations. How do I properly compute the actual residue of the integral above? Thanks for any suggestion.
EDIT
It was suggested by Robert Israel to rewrite
$$ \dfrac{1}{1+a+z_1/(z_1 + z_2)} = \dfrac{z_1 + z_2}{a z_1 + a z_2 + 2 z_1 + z_2}$$
and similar for the other term with $b$. This would give additional poles if these denominators vanish within $|z_1|=\epsilon$ and $|z_2|=\epsilon$. However, I meant these extra factors to be purely test functions that do not contribute any extra poles. Therefore, let us concentrate on the parameter region in $a$ and $b$ where the only poles in the above integrand are at $z_1=0$ and $z_2=0$. In this case the problem still stands unresolved.
(EDITED)
The poles at $z_1 = 0$ and $z_2=0$ are not the only ones inside the circles $|z|=\epsilon$. Note that $$ \dfrac{1}{z_1 z_2 (1+a+z_1/(z_1 + z_2)(1+b+z_2/(z_1+z_2))} = \dfrac{(z_1 + z_2)^2}{z_1 z_2 ((a+2) z_1 + (a+1) z_2)((b+1) z_1 + (b+2) z_2}$$ so (for fixed $z_2$ on the circle $|z|=\epsilon$) this has poles at $z_1 = 0$, $z_1 = -z_2 (a + 1)/(a + 2)$ and $z_1 = -z_2 (b+2)/(b+1)$. On the other hand, for fixed $z_1$ on the circle, it has poles at $z_2 = 0$, $z_2 = -z_1 (a+2)/(a+1)$ and $z_2 = -z_1 (b+1)/(b+2)$. Leaving out the cases $|a+2|=|a+1|$ and $|b+1|=|b+2|$ where there is a pole on the circle itself, one of $z_1 = -z_2 (a + 1)/(a + 2)$ and $z_2 = -z_1 (a+2)/(a+1)$ is inside the circle and the other is outside, and similarly for the poles involving $b$. The residues for poles inside the circle must be taken into account when doing the integrals.
For example, let's say $|a+1| > |a+2|$ and $|b+1| < |b+2|$. Then for the $z_1$ integral with fixed $z_2$ on the circle, we need only the pole at $z_1=0$:
$$ \eqalign{\frac{1}{2\pi i} \oint_{|z_1|=\epsilon} dz_1\; f(z_1, z_2) &= \text{res}(f(z_1,z_2); z_1 = 0) = \dfrac{1}{(b+2)(a+1) z_2}\cr \dfrac{1}{(2\pi i)^2} \oint_{|z_2|=\epsilon} dz_2 \oint_{|z_1|=\epsilon} dz_1\; f(z_1,z_2) &= \text{res}\left(\dfrac{1}{(b+2)(a+1)z_2}; z_2 = 0\right) = \dfrac{1}{(b+2)(a+1)}}$$ In the other order, for the $z_2$ integral with fixed $z_1$ on the circle we need the residues at $z_2 = 0$, $z_2 = -z_1 (a+2)/(a+1)$ and $z_2 = -z_1 (b+1)/(b+2)$.
$$ \eqalign{\frac{1}{2\pi i} \oint_{|z_2|=\epsilon} dz_2\; f(z_1, z_2) =& \text{res}(f(z_1,z_2); z_2 = 0) + \text{res}\left(f(z_1,z_2); z_2 = -z_1 \dfrac{a+2}{a+1}\right) \cr &+ \text{res}\left(f(z_1,z_2); z_2 = -z_1 \dfrac{b+1}{b+2}\right)\cr = {\frac {1}{z_{{1}} \left( b+1 \right) \left( a+2 \right) }}&+{\frac {1 }{z_{{1}} \left( a+2 \right) \left( a+1 \right) \left( a+b+3 \right) }}-{\frac {1}{z_{{1}} \left( b+2 \right) \left( b+1 \right) \left( a+b+3 \right) }}\cr &= \dfrac{1}{(b+2)(a+1) z_1}\cr \dfrac{1}{(2\pi i)^2} \oint_{|z_1|=\epsilon} dz_1 \oint_{|z_2|=\epsilon} dz_2\; f(z_1,z_2) &= \text{res}\left(\dfrac{1}{(b+2)(a+1) z_1}; z_1 = 0\right) = \dfrac{1}{(b+2)(a+1)} }$$ and there is no discrepancy.