We show that the entries of the stiffness matrix, associated with the $C^0$-piecewise linear finite element discretization of the hyper-singular integral fractional Laplacian (IFL) on rectangular meshes, can be simply expressed as one-dimensional integrals on a finite interval. Particularly, the FEM stiffness matrix on uniform meshes has a block-Toeplitz structure, so the matrix-vector multiplication can be implemented by FFT efficiently. The analytic integral representations not only allow for accurate evaluation of the entries, but also facilitate the study of some intrinsic properties of the stiffness matrix. For instance, we can obtain the asymptotic decay rate of the entries, so the “dense” stiffness matrix turns out to be “sparse” with an $\mathcal{O}(h^3)$ cutoff. We provide ample numerical examples of PDEs involving the IFL on rectangular or $L$-shaped domains to demonstrate the optimal convergence and efficiency of this semi-analytical approach. With this, we can also offer some benchmarks for the FEM on general meshes implemented by other means (e.g., for accuracy check and comparison when triangulation reduces to rectangular meshes).