  We present an algorithm for constructing the Wilson operator product
expansion (OPE) for perturbative interacting quantum field theory in general
Lorentzian curved spacetimes, to arbitrary orders in perturbation theory. The
remainder in this expansion is shown to go to zero at short distances in the
sense of expectation values in arbitrary Hadamard states. We also establish a
number of general properties of the OPE coefficients: (a) they only depend
(locally and covariantly) upon the spacetime metric and coupling constants, (b)
they satisfy an associativity property, (c) they satisfy a renormalization
group equation, (d) they satisfy a certain microlocal wave front set condition,
(e) they possess a ``scaling expansion''. The latter means that each OPE
coefficient can be written as a sum of terms, each of which is the product of a
curvature polynomial at a spacetime point, times a Lorentz invariant Minkowski
distribution in the tangent space of that point. The algorithm is illustrated
in an example.
