A general method for the calculation of the partition function of a canonical ensemble of noninteracting bound electrons is presented. It consists in a doubly recursive procedure with respect to the number of electrons and the number of orbitals. Contrary to existing approaches, this recursion relation contains no alternate summation of positive and negative numbers, which was the main source of numerical uncertainties. It is accompanied with a normalization of partition function through the determination of a free parameter consistent with the zeroth-order saddle-point approximation. The recursion relation allows one to calculate accurately partition functions for ions with a large number of orbitals, and is therefore important for calculations relying on the superconfiguration approximation.