  A new Monte Carlo method is proposed for fermion systems interacting with
classical degrees of freedom. To obtain a weight for each Monte Carlo sample
with a fixed configuration of classical variables, the moment expansion of the
density of states by Chebyshev polynomials is applied instead of the direct
diagonalization of the fermion Hamiltonian. This reduces a cpu time to scale as
$O(N_{\rm dim}^{2} \log N_{\rm dim})$ compared to $O(N_{\rm dim}^{3})$ for the
diagonalization in the conventional technique; $N_{\rm dim}$ is the dimension
of the Hamiltonian. Another advantage of this method is that parallel
computation with high efficiency is possible. These significantly save total
cpu times of Monte Carlo calculations because the calculation of a Monte Carlo
weight is the bottleneck part. The method is applied to the double-exchange
model as an example. The benchmark results show that it is possible to make a
systematic investigation using a system-size scaling even in three dimensions
within a realistic cpu timescale.
