Addition theorems for vector spherical harmonics require the determination of the Gaunt coefficients that appear in a linearization expansion of the product of two associated Legendre functions. This paper presents an algorithm for the efficient calculation of these coefficients through solving the most appropriate (lower triangular) linear system and derives all relevant recurrence relations needed in the calculation. This algorithm is also applicable to the calculation of the Clebsch-Gordan coefficients that are closely related to the Gaunt coefficients and are frequently encountered in the quantum theory of angular momentum. The new method described in this paper reduces the computing time to , compared to the existing formulation that is widely used. This new method can be applied to the calculation of both low- and high-degree coefficients, while the existing formulation works well only for low degrees.