  We describe a method for the numerical evaluation of normalized versions of
the associated Legendre functions $P_\nu^{-\mu}$ and $Q_\nu^{-\mu}$ of degrees
$0 \leq \nu \leq 1,000,000$ and orders $-\nu \leq \mu \leq \nu$ on the interval
$(-1,1)$. Our algorithm, which runs in time independent of $\nu$ and $\mu$, is
based on the fact that while the associated Legendre functions themselves are
extremely expensive to represent via polynomial expansions, the logarithms of
certain solutions of the differential equation defining them are not. We
exploit this by numerically precomputing the logarithms of carefully chosen
solutions of the associated Legendre differential equation and representing
them via piecewise trivariate Chebyshev expansions. These precomputed
expansions, which allow for the rapid evaluation of the associated Legendre
functions over a large swath of parameter domain mentioned above, are
supplemented with asymptotic and series expansions in order to cover it
entirely. The results of numerical experiments demonstrating the efficacy of
our approach are presented, and our code for evaluating the associated Legendre
functions is publicly available.
