  We present a method for studying the secular gravitational dynamics of
hierarchical multiple systems consisting of nested binaries, which is valid for
an arbitrary number of bodies and arbitrary hierarchical structure. We derive
the Hamiltonian of the system and expand it in terms of the -- assumed to be --
small ratios $x_i$ of binary separations. At the lowest nontrivial expansion
order (quadrupole order, second order in $x_i$), the Hamiltonian consists of
terms which, individually, depend on binary pairs. At higher orders, in
addition to terms depending on binary pairs, we also find terms which,
individually, depend on more than two binaries. In general, at order $n$ in
$x_i$, individual terms depend on at most $n-1$ binaries. We explicitly derive
the Hamiltonian including all terms up and including third order in $x_i$
(octupole order), and including the binary pairwise terms up and including
fifth order in $x_i$. These terms are orbit averaged, and we present a new
algorithm for efficiently solving the equations of motion. This algorithm is
highly suitable for studying the secular evolution of hierarchical systems with
complex hierarchies, making long-term integrations of such systems feasible. We
show that accurate results are obtained for multiplanet systems with semimajor
axis ratios as large as $\approx 0.4$, provided that high-order terms are
included. In addition to multiplanet systems with a single star, we apply our
results to multistar systems with multiple planets.
