  A high-order quadrature algorithm is presented for computing integrals over
curved surfaces and volumes whose geometry is implicitly defined by the level
sets of (one or more) multivariate polynomials. The algorithm recasts the
implicitly defined geometry as the graph of an implicitly defined, multi-valued
height function, and applies a dimension reduction approach needing only
one-dimensional quadrature. In particular, we explore the use of Gauss-Legendre
and tanh-sinh methods and demonstrate that the quadrature algorithm inherits
their high-order convergence rates. Under the action of $h$-refinement with $q$
fixed, the quadrature schemes yield an order of accuracy of $2q$, where $q$ is
the one-dimensional node count; numerical experiments demonstrate up to 22nd
order. Under the action of $q$-refinement with the geometry fixed, the
convergence is approximately exponential, i.e., doubling $q$ approximately
doubles the number of accurate digits of the computed integral. Complex
geometry is automatically handled by the algorithm, including, e.g.,
multi-component domains, tunnels, and junctions arising from multiple
polynomial level sets, as well as self-intersections, cusps, and other kinds of
singularities. A variety of numerical experiments demonstrates the quadrature
algorithm on two- and three-dimensional problems, including: randomly generated
geometry involving multiple high-curvature pieces; challenging examples
involving high degree singularities such as cusps; adaptation to simplex
constraint cells in addition to hyperrectangular constraint cells; and boolean
operations to compute integrals on overlapping domains.
