We discuss two algorithms for the numerical simulation of SU(3) lattice gauge theory with dynamical quarks. Both are based on the hybrid stochastic method of Duane and Kogut. They provide a relatively rapid evolution of the gauge fields through configuration space and good control of errors. One of the algorithms allows the simulation of arbitrary numbers of quarks. Tests of the algorithms are presented as well as initial data from a study of the thermodynamics of quarks and gluons with Kogut-Susskind fermions.