  A new stochastic number projection method is proposed. The component of the
BCS wave function corresponding to the right number of particles is obtained by
means of a Metropolis algorithm in which the weight functions are constructed
from the single-particle occupation probability. Either standard BCS or
Lipkin-Nogami probability distributions can be used, thus the method is
applicable for any pairing strength. The accuracy of the method is tested in
the computation of pairing energies of model and real systems.
