The production of jets in low Q 2 ep scattering (photoproduction) and in low Q 2 e + e − scattering (γγ scattering) allows for testing perturbative QCD and for measuring the proton and photon structure functions. This requires exact theoretical predictions for one-and two-jet cross sections. We describe the theoretical formalism, giving sufficient details, for calculating the direct and resolved processes in γp and γγ reactions in next-to-leading order QCD. We present the complete analytical results for the Born terms, the virtual, and the real corrections. To separate singular and regular regions of phase space we use the phase space slicing method with an invariant mass cut-off. In this way, all soft and collinear singularities are either canceled or absorbed into the structure functions. Using a flexible Monte Carlo program, we evaluate the cross sections numerically and perform various tests and comparisons with other calculations. We consider the scale dependence of our results and compare them to data from the experiments H1 and ZEUS at HERA and OPAL at LEP.