A hybrid UV/PBTG/SMCG method is developed based on the Numerical Solution of Maxwell Equation in 3D Simulations (NMM3D) of electromagnetic wave scattering by Gaussian dielectric random rough surfaces with exponential correlation functions. In 3D simulations, the height function z=f(x,y) of the rough surfaces vary in both two horizontal directions. The method inherits the advantages of the UV method, the Physics Based Two Grid (PBTG) method, and the Sparse Matrix Canonical Grid (SMCG) method. The method is implemented in parallel computing. Numerical simulations show that the hybrid method is several times more efficient than the previous approaches. With the improved computational efficiency, we are able to compute cases with surface areas up to 32 by 32 square wavelengths and rms heights up to half wavelengths. To verify the accuracies of the results, we did 5 tests: convergence with discrete samplings, convergence with realizations, convergence with surface sizes, energy conservation for each realization, and reciprocity for each realization. New physical results including backscattering enhancements were obtained.