This paper describes an efficient algorithm for computing steady two-dimensional surface gravity waves in irrotational motion. The algorithm complexity is $O(N\log N)$, $N$ being the number of Fourier modes. This feature allows the arbitrary precision computation of waves in arbitrary depth, i.e. it works efficiently for Stokes, cnoidal and solitary waves, even for quite large steepnesses, up to approximately 99 % of the maximum steepness for all wavelengths. In particular, the possibility to compute very long (cnoidal) waves accurately is a feature not shared by other algorithms and asymptotic expansions. The method is based on conformal mapping, the Babenko equation rewritten in a suitable way, the pseudo-spectral method and Petviashvili iterations. The efficiency of the algorithm is illustrated via some relevant numerical examples. The code is open source, so interested readers can easily check the claims, use and modify the algorithm.