In this work we revisit the problem of the dynamical stability of hierarchical triple systems with applications to circumbinary planetary orbits. We derive critical semimajor axes based on simulating and analyzing the dynamical behavior of 3 × 108 binary star–planet configurations. For the first time, three-dimensional and eccentric planetary orbits are considered. We explore systems with a variety of binary and planetary mass ratios, binary and planetary eccentricities from 0 to 0.9, and orbital mutual inclinations ranging from 0° to 180°. Planetary masses range between the size of Mercury and the lower fusion boundary (approximately 13 Jupiter masses). The stability of each system is monitored over 106 planetary orbital periods. We provide empirical expressions in the form of multidimensional, parameterized fits for two borders that separate dynamically stable, unstable, and mixed zones. In addition, we offer a machine learning model trained on our data set as an alternative tool for predicting the stability of circumbinary planets. Both the empirical fits and the machine learning model are tested for their predictive capabilities against randomly generated circumbinary systems with very good results. The empirical formulae are also applied to the Kepler and TESS circumbinary systems, confirming that many planets orbit their host stars close to the stability limit of those systems. Finally, we present a REST application programming interface with a web-based application for convenient access to our simulation data set.