  An improved algorithm is proposed for Monte Carlo methods to study fermion
systems interacting with adiabatical fields. To obtain a weight for each Monte
Carlo sample with a fixed configuration of adiabatical fields, a series
expansion using Chebyshev polynomials is applied. By introducing truncations of
matrix operations in a systematic and controlled way, it is shown that the cpu
time is reduced from O(N^3) to O(N) where N is the system size. Benchmark
results show that the implementation of the algorithm makes it possible to
perform systematic investigations of critical phenomena using system-size
scalings even for an electronic model in three dimensions, within a realistic
cpu timescale.
