In the spherical harmonic model for the gravitational potential and its derivatives, singularities are encountered at the poles. These singularities, caused by the use of spherical coordinates, can be removed with direction cosines and “derived” Legendre functions. Using the definition and properties of derived Legendre functions, a set of recursion formulas may be formed for their computation. However, by imposing a set of algorithmic constraints, the set of suitable recursion formulas is reduced to six. Finally, a numerical stability requirement reduces the number of suitable relationships from six to two. The computational precision associated with the acceptable recursions has been examined for degree and order up to 180, as required for some forthcoming satellite missions. The precision of the gravitation force computation obtained with the acceptable recursions was demonstrated to be near the limiting precision of the machine.