We present the publicly available program NGluon allowing the numerical evaluation of primitive amplitudes at one-loop order in massless QCD. The program allows the computation of one-loop amplitudes for an arbitrary number of gluons. The focus of the present article is the extension to one-loop amplitudes including an arbitrary number of massless quark pairs. We discuss in detail the algorithmic differences to the pure gluonic case and present cross checks to validate our implementation. The numerical accuracy is investigated in detail.