Computing electromagnetic properties of large-scale finite periodic structures (LFPSs) via the multilevel fast multipole algorithm (MLFMA) often suffers from heavy time and memory consumption. The main reason is that complex structures and high-permittivity materials in LFPSs need dense meshes, leading to a prohibitively dense near-field interaction (NFI) matrix. To overcome this drawback, during the MLFMA initialization, each leaf box is enforced to exactly encapsulate an identical periodic structure of the LFPSs with consistent mesh topology and basis function definitions. The NFI matrix can then be divided into several repeatable submatrices, while only those unique ones require computation and storage, thereby reducing the filling time and memory consumption of the NFI matrix by severalfold without sacrificing accuracy. Additionally, when the leaf boxes are overcrowded with basis and test functions, the matrix-vector products during the iterative solution are enhanced through dual accelerations: interpolative decomposition (ID) compresses the low-rank NFI matrix, while spherical harmonic expansion (SE) replaces numerical quadrature in k-space with spherical harmonic summations for far-field interactions. Two LFPSs with commonly rectangular and triangular lattice arrangements are computed to verify the efficiency and accuracy of the proposed method.