The two-loop QCD corrections to vector boson pair production at hadron colliders involve a new class of Feynman integrals: two-loop four-point functions with two off-shell external legs. We describe their reduction to a small set of master integrals by solving linear relations among them. We then use differential equations in the external invariants to compute all master integrals that are relevant to planar Feynman amplitudes. Our results are expressed analytically in terms of generalized harmonic polylogarithms. The calculation relies heavily on techniques that exploit the algebraic structure of these functions, which we describe in detail.