An accelerated Integral Equations (IE) field solver for determining scattered fields from electrically large electromagnetic metasurfaces utilizing Fast Multipole Method (FMM) is proposed and demonstrated in 2D. In the proposed method, practical general metasurfaces are expressed using an equivalent zero thickness sheet model described using surface susceptibilities, and where the total fields around it satisfy the Generalized Sheet Transition Conditions (GSTCs). While the standard IE-GSTC offers fast field computation compared to other numerical methods, it is still computationally demanding when solving electrically large problems, with a large number of unknowns. Here we accelerate the IE-GSTC method using the FMM technique by dividing the current elements on the metasurface into near- and far-groups, where either the rigorous or approximated Green’s function is used, respectively, to reduce the computation time without losing solution accuracy. Using numerical examples, the speed improvement of the FMM IE-GSTC method O(N^1.5) over the standard IE-GSTC O(N³) method is confirmed. Finally, the usefulness of the FMM IE-GSTC is demonstrated by applying it to solve electromagnetic propagation inside an electrically large radio environment with strategically placed metasurfaces to improve signal coverage in blind areas, where a standard IE-GSTC solver would require prohibitively large computational resources and long simulation times.