We describe a surface integral-equation (SIE) method suitable for computation of electromagnetic fields scattered by 2D-periodic high-permittivity and plasmonic scatterers. The method makes use of fast evaluation of the 2D-quasi-periodic Green function (2D-QPGF) and its gradient using a tabulation technique in combination with tri-linear interpolation. In particular we present a very efficient technique to create the look-up tables for the 2D-QPGF and its gradient where we use to our advantage that it is very effective to simultaneously compute the QPGF and its gradient, and to simultaneously compute these values for the case in which the role of source and observation point are interchanged. We use the Ewald representation of the 2D-QPGF and its gradient to construct the tables with pre-computed values. Usually the expressions for the Ewald representation of the 2D-QPGF and its gradient are presented in terms of the complex complementary error function but here we give the expressions in terms of the Faddeeva function enabling efficient use of the dedicated algorithms to compute the Faddeeva function. Expressions are given for both lossy and lossless medium parameters and it is shown that the expression for the lossless case can be evaluated twice as fast as the expression for the lossy case. Two case studies are presented to validate the proposed method and to show that the time required for computing the method of moments (MoM) integrals that require evaluation of the 2D-QPGF becomes comparable to the time required for computing the MoM integrals that require evaluation of the aperiodic Green function.