  We present a polynomial hybrid Monte Carlo (PHMC) algorithm for lattice QCD
with odd numbers of flavors of O(a)-improved Wilson quark action. The algorithm
makes use of the non-Hermitian Chebyshev polynomial to approximate the inverse
square root of the fermion matrix required for an odd number of flavors. The
systematic error from the polynomial approximation is removed by a noisy
Metropolis test for which a new method is developed. Investigating the property
of our PHMC algorithm in the N_f=2 QCD case, we find that it is as efficient as
the conventional HMC algorithm for a moderately large lattice size (16^3 times
48) with intermediate quark masses (m_{PS}/m_V ~ 0.7-0.8). We test our
odd-flavor algorithm through extensive simulations of two-flavor QCD treated as
an N_f = 1+1 system, and comparing the results with those of the established
algorithms for N_f=2 QCD. These tests establish that our PHMC algorithm works
on a moderately large lattice size with intermediate quark masses (16^3 times
48, m_{PS}/m_V ~ 0.7-0.8). Finally we experiment with the (2+1)-flavor QCD
simulation on small lattices (4^3 times 8 and 8^3 times 16), and confirm the
agreement of our results with those obtained with the R algorithm and
extrapolated to a zero molecular dynamics step size.
