Contact of rough surfaces is of prime importance in the study of friction and wear. Numerical simulations are well suited for this non-linear problem, but natural surfaces being fractal [1], they have high discretization requirements. There is therefore a need for efficient numerical methods able to deal with non-linearities. With finite elements [2], it is necessary to discretize the volume of the solids in contact, increasing the computational cost. Efficient optimization methods [3, 4] for the contact problem with a boundary-element formulation based on the FFT algorithm [5] allow an economy of computational resources, as only the boundary of the solids is discretized. Extension of a boundary-element method for elasto-plastic contact is feasible [6], on the basis of a boundary-domain integral formulation. We present here a new version of this formulation that takes advantage of the FFT for a reduction of both algorithmic and memory complexity. The method is based on our new derivation of the Mindlin [7] fundamental solution using a partial Fourier transform, similarly to [8]. This allows the computation of displacement gradients necessary for applying a Newton-Raphson iterative solver with a consistent tangent operator to the boundary-domain integral formulation [9]. We present an adaptation of this treatment that takes advantage of the FFT. The solver gives the plastic deformations, which are used to compute a new contact solution [6]. With this new method, we compare to models of elastic rough contact [10, 11] and elastic-saturated rough contact [12, 13], and quantify the effects of elasto-plasticity on the contact area and the distribution of contact clusters, the latter playing an important role in the characterization of wear of surfaces [14].