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. There are a number of algorithms for these functions published since 1960 but none of them satisfy our requirements. 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 10E-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.