Accurate simulations of molecular quantum dynamics are crucial for understanding numerous natural processes and experimental results. Yet, such high-accuracy simulations are challenging even for relatively simple systems where the Born-Oppenheimer approximation is valid. Worse still, due to the ubiquity of conical intersections in molecules, it is often necessary to go beyond this celebrated approximation and employ even more sophisticated methods that can describe nonadiabatic processes involving multiple coupled electronic states. Geometric integrators provide a way to simulate nonadiabatic quantum dynamics with high accuracy while exactly conserving geometric properties, such as norm, energy, symplecticity, and time reversibility. However, the exact conservation of these properties typically leads to a higher computational cost when compared with non-geometric integrators of the same order of accuracy. We remedy this lack of efficiency by employing various composition methods to obtain high-order geometric integrators, competitive in efficiency even with some non-geometric integrators. To overcome the limited applicability of the popular split-operator algorithms, we present geometric integrators that can solve the time-dependent Schroedinger equation efficiently regardless of the form of the Hamiltonian. Employing these integrators, we systematically compare the different representations of the molecular Hamiltonian for their suitability for simulating nonadiabatic dynamics at a conical intersection and find that the rarely used exact quasidiabatic Hamiltonian, which includes the often-neglected residual nonadiabatic couplings, yields the most accurate results. Accordingly, we present a method for quantifying the validity of ignoring these residual couplings and show that depending on the system, initial state, and employed quasidiabatization scheme, neglecting the residual couplings can indeed result in inaccurate dynamics. The computational cost of exact quantum simulations quickly becomes prohibitively high as the system size increases. For high-dimensional simulations, one can either optimize the available basis set, e.g., by employing adaptive phase space grids or, alternatively, employ a more approximate approach, such as mixed quantum-classical Ehrenfest dynamics, that provides a useful qualitative picture. Although Ehrenfest dynamics and exact quantum dynamics simulated on an adaptive grid may appear unrelated, we show that, surprisingly, the high-order geometric integrators for these two problems can, in fact, be obtained by using the same splitting and composition methods.