We give a detailed analysis of the origin of spurious divergences and finite steps that have been recently identified in particle-number restoration calculations within the nuclear energy density functional framework. We isolate two distinct levels of spurious contributions to the energy. The first one is encoded in the definition of the basic energy density functional itself whereas the second one relates to the canonical procedure followed to extend the use of the energy density functional to multi-reference calculations. The first level of spuriosity relates to the long-known self-interaction problem and to the newly discussed self-pairing interaction process which might appear when describing paired systems with energy functional methods using auxiliary reference states of Bogoliubov or BCS type. A minimal correction to the second level of spuriosity to the multi-reference nuclear energy density functional proposed in [D. Lacroix, T. Duguet, M. Bender, arXiv:0809.2041] is shown to remove completely the anomalies encountered in particle-number restored calculations. In particular, it restores sum-rules over (positive) particle numbers that are to be fulfilled by the particle-number-restored formalism. The correction is found to be on the order of several hundreds of keVs up to about 1 MeV in realistic calculations, which is small compared to the total binding energy, but often accounts for a substantial percentage of the energy gain from particle-number restoration and is on the same energy scale as the excitations one addresses with multi-reference energy density functional methods.