Conventional estimators of the anisotropic power spectrum and two-point correlation function (2PCF) adopt the `Yamamoto approximation', fixing the line-of-sight of a pair of galaxies to that of just one of its members. Whilst this is accurate only to first-order in the characteristic opening angle $\theta_\max$, it allows for efficient implementation via Fast Fourier Transforms (FFTs). This work presents practical algorithms for computing the power spectrum and 2PCF multipoles using pairwise lines-of-sight, adopting either the galaxy midpoint or angle bisector definitions. Using newly derived infinite series expansions for spherical harmonics and Legendre polynomials, we construct estimators accurate to arbitrary order in $\theta_\max$, though note that the midpoint and bisector formalisms themselves differ at fourth order. Each estimator can be straightforwardly implemented using FFTs, requiring only modest additional computational cost relative to the Yamamoto approximation. We demonstrate the algorithms by applying them to a set of realistic mock galaxy catalogs, and find both procedures produce comparable results for the 2PCF, with a slight preference for the bisector power spectrum algorithm, albeit at the cost of greater memory usage. Such estimators provide a useful method to reduce wide-angle systematics for future surveys.