Multivariate Residue Theorem?

There are several problems with residues of functions of several complex variables, but they are not insurmountable.

Let's begin by thinking about the one-variable case. To get something that can be generalized to several variables, we have to think about the one-variables residues in a particular way. First of all, we should not think of residues as an operation attaching numbers to functions, but as an operation on differential forms.

In one complex variable, we normally integrate around curves, most often curves that bound a domain in $\mathbb{C}$, and the proper thing to integrate along a curve is a $1$-form. Differential $1$-forms come in two flavors: $(1,0)$-forms, i.e. things that locally can be written $f(z,\bar z)\,dz$ and $(0,1)$-forms of the type $f(z,\bar z)\,d\bar z$. In a first course on complex analysis, we usually only find the first version, and most often when $f$ is holomorphic, i.e. when ''$f$ only depends on $z$, not on $\bar z$''.

It turns out that a good way to think of residues is as something that associates a new "form", (technically a current) to a $1$-form $f(z)\,dz$. This new form, let us call it a residue current should be supported on the singularity set of $f$, and we want to have a representation formula

$$\int_\gamma f(z)\,dz = \int_X \operatorname{res}(f).$$

In the simplest case, when $f$ has a single simple pole inside $\gamma$, the residue of $f$ should be a point mass located at the singularity; more precisely $\operatorname{res}(\frac{f(z)}{z-p}) = 2\pi i f(p) \delta_p$, where $p$ is the single pole. With this interpretation, we get the usual Cauchy's integral formula:

$$\int_\gamma \frac{f(z)}{z-p}\,dz = \int_{p} \operatorname{res}(\frac{f(z)}{z-p}) = \int_{p} 2\pi i f(p)\,d\delta_p(z) = 2\pi i f(p) .$$

In several variables, things get more complicated. Curves do not bound domains, and singularities of meromorphic functions are not compact, and not $0$-dimensional. If we start with a $(n,n-p)$-form, we want to associate a current, supported on the singularity of the form in such a way that we get an integral representation form similar to the one above.

For forms having singularities of particularly simple type (keyword: complete intersections), it's been known for a long time how to do this. Coleff and Herrera wrote an influential paper in 1978 about this. For more complicated singularities, there have been a number of improvements and generalizations, due to (among others) Passare, Tsikh, Andersson, Samuelsson, Wulcan just to mention a few. It is still an active research topic to find the "best" version of residue currents, and it requires tools such from algebraic geometry (including Hironaka's resolution of singularities) as well as commutative algebra (various types of ideal decompositions).

You can get a higher-dimensional version of the residue theorem by applying the one-dimensional version repeatedly. For example, let $f(x_1,x_2)$ be analytic on $\bar\Omega_1 \times \bar \Omega_2$, and let $x_1 \in \Omega_1$, $x_2 \in \Omega_2$. Then, $$ -\frac{1}{4\pi^2} \int_{\partial\Omega_1} \int_{\partial\Omega_2} \frac{f(t_1,t_2)}{(t_1 - x_1) \, (t_2 - x_2)} \, dt_1 \, dt_2 = \frac{1}{2\pi i} \int_{\partial\Omega_1} \frac{f(t_1,x_2)}{t_1 - x_1} \, dt_1 = f(x_1,x_2) . $$ The situation can get messy, however, if your pole does not have tensor-product form. For example, if we replace $(t_1 - x_1) \, (t_2 - x_2)$ with $t_1 - t_2 + x$, we get $$ -\frac{1}{4\pi^2} \int_{\partial\Omega_1} \int_{\partial\Omega_2} \frac{f(t_1,t_2)}{t_1 - t_2 + x} \, dt_1 \, dt_2 = - \frac{1}{2\pi i} \int_{\partial\Omega_1} f(t_1,t_1 + x) \, dt_1 , $$ and now you have to think about residues of $g(x_1) := f(x_1,x_1+x)$ to determine the last integral.