Explicitly-correlated F12 methods are becoming the first choice for high-accuracy molecular orbital calculations, and can often achieve chemical accuracy with relatively small gaussian basis sets. In most calculations, the many three-and four-electron integrals that formally appear in the theory are avoided through judicious use of resolutions of the identity (RI). However, in order not to jeopardize the intrinsic accuracy of the F12 wave function, the associated RI auxiliary basis set must be large. Here, inspired by the Head-Gordon-Pople (HGP) and PRISM algorithms for two-electron integrals, we present an algorithm to compute directly three-electron integrals over gaussian basis functions and a very general class of three-electron operators, without invoking RI approximations. A general methodology to derive vertical, transfer and horizontal recurrence relations is also presented.