We present a polynomial Hybrid Monte Carlo (PHMC) algorithm as an exact simulation algorithm with dynamical Kogut-Susskind fermions. The algorithm uses a Hermitian polynomial approximation for the fractional power of the KS fermion matrix. The systematic error from the polynomial approximation is removed by the Kennedy-Kuti noisy Metropolis test so that the algorithm becomes exact at a finite molecular dynamics step size. We performed numerical tests withNf=2 case on several lattice sizes. We found that the PHMC algorithm works on a moderately large lattice of 16 4 atβ=5.7, m=0.02 (mps/mv≈0.69) with a reasonable computational time.