  We present an algorithmic approach to the verification of identities on
multiple theta functions in the form of products of theta functions
$[(-1)^{\delta}a_1^{\alpha_1}a_2^{\alpha_2}\cdots a_r^{\alpha_r}q^{s};
q^{t}]_\infty$, where $\alpha_i$ are integers, $\delta=0$ or $1$, $s\in
\mathbb{Q}$, $t\in \mathbb{Q}^{+}$, and the exponent vectors
$(\alpha_1,\alpha_2,\ldots,\alpha_r)$ are linearly independent over
$\mathbb{Q}$. For an identity on such multiple theta functions, we provide an
algorithmic approach for computing a system of contiguous relations satisfied
by all the involved multiple theta functions. Using Stanley's Lemma on the
fundamental parallelepiped, we show that a multiple theta function can be
determined by a finite number of its coefficients. Thus such an identity can be
reduced to a finite number of simpler relations. Many classical multiple theta
function identities fall into this framework, including Riemann's addition
formula and the extended Riemann identity.
