Associated Legendre polynomials and spherical harmonics are central to calculations in many fields of science and mathematics – not only chemistry but computer graphics, magnetic, seismology and geodesy. In this paper, we present a comprehensive review of algorithms in the literature and, based on them, propose an efficient and accurate code for quantum chemistry. Our requirements are to efficiently calculate these functions for all non-negative integer degrees and orders up to a given number (<= 1000) and the absolute or the relative error of each calculated value should not exceed 10^-10. We achieve this by normalizing the polynomials, employing efficient and stable recurrence relations, and precomputing coefficients. The algorithm presented here is straightforward and may be used in other areas of science.