  The production of jets in low $Q^2$ $ep$ scattering (photoproduction) and in
low $Q^2$ $e^+e^-$ scattering ($\gamma\gamma$ 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 $\gamma p$ and $\gamma\gamma$
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.
