  We describe and implement an exact, flexible, and computationally efficient
algorithm for joint component separation and CMB power spectrum estimation,
building on a Gibbs sampling framework. Two essential new features are 1)
conditional sampling of foreground spectral parameters, and 2) joint sampling
of all amplitude-type degrees of freedom (e.g., CMB, foreground pixel
amplitudes, and global template amplitudes) given spectral parameters. Given a
parametric model of the foreground signals, we estimate efficiently and
accurately the exact joint foreground-CMB posterior distribution, and therefore
all marginal distributions such as the CMB power spectrum or foreground
spectral index posteriors. The main limitation of the current implementation is
the requirement of identical beam responses at all frequencies, which restricts
the analysis to the lowest resolution of a given experiment. We outline a
future generalization to multi-resolution observations. To verify the method,
we analyse simple models and compare the results to analytical predictions. We
then analyze a realistic simulation with properties similar to the 3-yr WMAP
data, downgraded to a common resolution of 3 degree FWHM. The results from the
actual 3-yr WMAP temperature analysis are presented in a companion Letter.
