Gravity data observed on or reduced to the ellipsoid are preferably represented using ellipsoidal harmonics instead of spherical harmonics. Ellipsoidal harmonics, however, are difficult to use in practice because the computation of the associated Legendre functions of the second kind that occur in the ellipsoidal harmonic expansions is not straightforward. Jekeli’s renormalization simplifies the computation of the associated Legendre functions. We extended the direct computation of these functions—as well as that of their ratio—up to the second derivatives and minimized the number of required recurrences by a suitable hypergeometric transformation. Compared with the original Jekeli’s renormalization the associated Legendre differential equation is fulfilled up to much higher degrees and orders for our optimized recurrences. The derived functions were tested by comparing functionals of the gravitational potential computed with both ellipsoidal and spherical harmonic syntheses. As an input, the high resolution global gravity field model EGM2008 was used. The relative agreement we found between the results of ellipsoidal and spherical syntheses is 10⁻¹⁴, 10⁻¹² and 10⁻⁸ for the potential and its first and second derivatives, respectively. Using the original renormalization, this agreement is 10⁻¹², 10⁻⁸ and 10⁻⁵, respectively. In addition, our optimized recurrences require less computation time as the number of required terms for the hypergeometric functions is less.