Reconstructions of past solar activity based on cosmogenic radioisotopes have reavealed that the Sun spends a significant fraction (≈20 %) of its time in aperiodically recurring states of so-called Grand Minima or Grand Maxima, namely epochs of strongly supressed and markedly above-average levels of magnetic activity, respectively. The physical origin of these episodes is not yet understood. In this article we present a dual-dynamo model of the solar cycle, combining a dominant dynamo based on differential-rotation shear and surface decay of bipolar active regions, and a weak, deep-seated turbulent dynamo. The resulting dynamo simulations are found to exhibit the equivalent of observed Grand Minima and Maxima. By adjusting the magnitude and saturation level of the secondary turbulent dynamo, we can reproduce well the duration and waiting-time distributions of Grand Minima and Maxima inferred from the cosmogenic-isotope record. The exit from Grand Minima episodes is typically characterized by strong hemispheric asymmetries, in agreement with sunspot observations during the 1645 - 1715 Maunder Minimum. In these simulations, Grand Maxima can be unambiguously identified as a distinct dual-dynamo state resulting from constructive interference between the two dynamos mechanisms operating within the simulation. This interaction leads to the autonomous production of long quasi-periodicities in the millennial range, commensurate with the Halstatt cycle. Such a quasi-periodic modulation, readily produced through dynamical backreaction on large-scale flows in non-kinematic dynamo models, is quite uncommon in a purely kinematic solar-cycle model such as the one developed herein. We argue that these long periodicities are set by the long diffusion time of magnetic field accumulating in the stable layers underlying the turbulent convection zone.