We propose an exact simulation algorithm for lattice QCD with dynamical Kogut–Susskind fermion in which the Nf-flavor fermion operator is defined as the Nf/4th root of the Kogut–Susskind (KS) fermion operator. The algorithm is an extension of the Polynomial Hybrid Monte Carlo (PHMC) algorithm to KS fermions. The fractional power of the KS fermion operator is approximated with a Hermitian Chebyshev polynomial, with which we can construct an algorithm for any number of flavors. The error which arises from the approximation is corrected by the Kennedy–Kuti noisy Metropolis test. Numerical simulations are performed for the two-flavor case for several lattice parameters in order to confirm the validity and the practical feasibility of the algorithm. In particular tests on a 16 4 lattice with a quark mass corresponding to mPS/mVnot, vert, similar0.68 are successfully accomplished. We conclude that our algorithm provides an attractive exact method for dynamical QCD simulations with KS fermions.