  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.
