We present a new high-order numerical method for the solution of high-frequency scattering problems from rough surfaces in three dimensions. The method is based on the rigorous solution of appropriate integral equations in the high-frequency regime, in a manner that bypasses the need to resolve the fields on the scale of the wavelength of radiation. Indeed, inspired by prior work in two dimensions [Bruno, Sei & Caponi, Radio Sci. 37 (2002)], we seek a solution in the form of a slow modulation of the incoming radiation, and we choose a series expansion in inverse powers of the wavenumber to represent the unknown slowly varying envelope. As we show, this framework can be made to yield an efficiently computable recursion for the terms in the series, which is similar in nature to its two dimensional counterpart. The actual derivation of the recursion, however, differs from that of the two-dimensional case as the stronger singularity of the three-dimensional free-space Green's function demands a more careful treatment. Our approach then is based on an "analytic" resolution of this singularity (effected via a change of independent variables), and on a subsequent explicit asymptotic treatment of the resulting oscillatory integrals, which allows us to derive (convergent) expansions of arbitrarily high order in wavelength. As we demonstrate, the resulting algorithms generally provide a very significant improvement over classical (e.g. Kirchhoff's) approximations in both accuracy and applicability and they can, in fact, efficiently produce results with full double-precision accuracy for configurations of   practical interest.