Quasi-2D Coulomb systems are of fundamental importance and have attracted much attention in many areas nowadays. Their reduced symmetry gives rise to interesting collective behaviors, but also brings great challenges for particle-based simulations. We propose a novel algorithm framework to address the O(N^2) simulation complexity associated with the long-range nature of Coulomb interactions. First, we introduce an efficient Sum-of-Exponentials (SOE) approximation for the long-range kernel associated with Ewald splitting, achieving uniform convergence in terms of inter-particle distance, which reduces the complexity to O(N^7/5). We then introduce a random batch sampling method in the periodic dimensions, the stochastic approximation is proven to be both unbiased and with reduced variance via a tailored importance sampling strategy, further reducing the computational cost to O(N). The performance of our algorithm is demonstrated via varies numerical examples. Notably, it achieves a speedup of 2 ∼3 orders of magnitude comparing with Ewald2D method, enabling molecular dynamics (MD) simulations with up to 10^6 particles on a single core. The present approach is therefore well-suited for large-scale particle-based simulations of Coulomb systems under confinement, making it possible to investigate the role of Coulomb interaction in many practical situations.