Suppose I have two discrete-support random variables, $X$ and $Y$.
They have joint CDF $F(X,Y)$.
If I want to find
$\Pr(a \leq X \leq b , c \leq Y \leq d)$.
It is obviously not:
$F(b ,d)-F(a-1 ,c-1)$.
What is it then, in simplest terms with respect to the CDF?
I believe, for two dimensions, it is:
$F(b ,d)-F(a-1 ,d)-F(b,c-1)+F(a-1,c-1)$.
What is it for an arbitrary $d$ dimensions?
The distribution function (CDF) of a random variable $X=(X_1,X_2,\ldots, X_n)$ with values in $\mathbb{R}^n$ is defined to be
$$F_X(x_1,x_2,\ldots, x_n) = \Pr(X_1\le x_1,\, X_2\le x_2,\ \ldots,\ X_n\le x_n).$$
To find the chance that $X$ lies within some nonempty (half-open) parallelipiped $$\mathcal{I} = (i_1,j_1] \times (i_2,j_2] \times \cdots \times (i_n,j_n]$$
(or equivalently, to find the chance that simultaneously $i_k \lt X_k \le j_k$ for $k=1,2,\ldots,n$), simplify the problem by choosing coordinate axes centered on $(i_1,i_2,\ldots,i_n)$ with unit distances $j_k-i_k.$ This makes all the $i_k$ equal to $0$ and all the $j_k$ equal to $1:$ $\mathcal I$ becomes the unit cube in $\mathbb{R}^n.$ (This is really just a notational convenience rather than a change of any mathematical substance.)
By partitioning the last coordinate into non-positive values and values less than or equal to $1,$ observe that the unit cube -- merely as a set -- can be expressed as the set difference
$$\mathcal{I} = \mathcal{I}^{(1)}\setminus\mathcal{I}^{(0)}$$
where (for any real number $s$) $$\mathcal{I}^{(s)} = (0,1]\times(0,1]\times \cdots \times (0,1] \times (-\infty, s].$$
The axioms of probability (specifically, that probabilities of non-intersection regions will add) imply
$$\eqalign{ \Pr(X\in\mathcal{I}) &= \Pr(X\in\mathcal{I}^{(1)}) - \Pr(X\in\mathcal{I}^{(0)}) \\ &= F(1,1,\ldots,1,1) - F(1,1,\ldots,1,0). }$$
Repeating this reduction at coordinates $n-1,$ $n-2, ...$ down to $1$ will double the number of terms at each iteration. The resulting combination therefore has $2^n$ terms. All terms are of the form $\pm F(s_1,s_2,\ldots,s_n)$ with $s_j\in\{0,1\},$ and therefore are distinct: there is no algebraic simplification. The sign is determined by the number of zeros in the arguments of $F$: negative for an odd number, positive for an even number.
This result can be written as an iterated sum over $n$ variables,
$$\Pr(X\in (0,1]^n) = \sum_{s_1\in\{0,1\}}\sum_{s_2\in\{0,1\}}\cdots \sum_{s_n\in\{0,1\}}(-1)^{n-\sum_{j=1}^n s_j}\ F(s_1,s_2, \ldots, s_n).$$
In a more succinct and memorable vector notation, where $\mathbf{s} = (s_1,s_2,\ldots,s_n)$ and $|\mathbf{s}| = s_1+s_2+\cdots+s_n,$ we may write
$$\Pr(X\in (0,1]^n) = \sum_{\mathbf{s}\in\{0,1\}^n} (-1)^{n-|\mathbf{s}|} \ F(\mathbf{s}).$$
When $n=3,$ for instance, this is
$$\Pr(X\in \mathcal (0,1]^3) = F(1,1,1)-F(1,1,0)-F(1,0,1)-F(0,1,1)+F(1,0,0)+F(0,1,0)+F(0,0,1)-F(0,0,0).$$
In terms of the original coordinates this would be written
$$\Pr(X\in \mathcal I) = F(j_1,j_2,j_3)-F(j_1,j_2,i_3)-F(j_1,i_2,j_3)-F(i_1,j_2,j_3)+F(j_1,i_2,i_3)+F(i_1,j_2,i_3)+F(i_1,i_2,j_3)-F(i_1,i_2,i_3).$$
Finally, when $X$ is supported on the integral lattice $\mathbb{Z}^n$ (of vectors all of whose coordinates are integers), because
$$\Pr(X_j \lt x) = \Pr(X_j \le x-1),$$
you may replace each $i_k$ by $i_k-1$ to find the chance that $X$ lies in the closed parallelipiped $$\bar{\mathcal{I}} = [i_1,j_1]\times \cdots \times [i_n,j_n].$$
