A high-order quadrature rule is constructed for the evaluation of Laplace single and double layer potentials and their normal derivatives on smooth surfaces in three dimensions. The construction starts from a harmonic approximation of the density on each patch that admits a natural harmonic polynomial extension in a volumetric neighborhood of the patch. Then by the general Stokes theorem, singular and nearly singular surface integrals are reduced to line integrals preserving the singularity of the kernel. These singularity-preserving line integrals can be evaluated semi-analytically by the singularity-swap quadrature involving a complex-valued parameter root finding for the given target and polynomial approximation of the smooth part on the boundary of the patch. In other words, the evaluation of singular and nearly singular surface integrals is reduced to function evaluations on the vertices on the boundary of each patch. The recursive reduction quadrature largely removes adaptive integration that is needed in most existing high-order quadratures for singular and nearly singular surface integrals, leading to improved efficiency and robustness.