  The mean values of a many-body Hamiltonian including a proton-neutron pairing
term and matrix elements of one-, two- and four-body operators within a basis
of particle number projected BCS states, are analytically expressed in terms of
a single function Q(N) depending on the number of particles, $N$. The function
Q(N) is calculated using a recursion in $N$ in which the shells and the BCS
angles are kept the same for any step of iteration. An illustrative example is
numerically considered in a restricted single particle space. Some specific
features for the standard BCS, the projection after variation approach as well
as for the variation after projection formalism, are pointed out.
