<p style='text-indent:20px;'>In this paper, we propose a new class of operator factorization methods to discretize the integral fractional Laplacian <inline-formula><tex-math id="M1">\begin{document}$ (- \Delta)^\frac{{ \alpha}}{{2}} $\end{document}</tex-math></inline-formula> for <inline-formula><tex-math id="M2">\begin{document}$ \alpha \in (0, 2) $\end{document}</tex-math></inline-formula>. One main advantage is that our method can easily increase numerical accuracy by using high-degree Lagrange basis functions, but remain its scheme structure and computer implementation unchanged. Moreover, it results in a symmetric (multilevel) Toeplitz differentiation matrix, enabling efficient computation via the fast Fourier transforms. If constant or linear basis functions are used, our method has an accuracy of <inline-formula><tex-math id="M3">\begin{document}$ {\mathcal O}(h^2) $\end{document}</tex-math></inline-formula>, while <inline-formula><tex-math id="M4">\begin{document}$ {\mathcal O}(h^4) $\end{document}</tex-math></inline-formula> for quadratic basis functions with <inline-formula><tex-math id="M5">\begin{document}$ h $\end{document}</tex-math></inline-formula> a small mesh size. This accuracy can be achieved for any <inline-formula><tex-math id="M6">\begin{document}$ \alpha \in (0, 2) $\end{document}</tex-math></inline-formula> and can be further increased if higher-degree basis functions are chosen. Numerical experiments are provided to approximate the fractional Laplacian and solve the fractional Poisson problems. It shows that if the solution of fractional Poisson problem satisfies <inline-formula><tex-math id="M7">\begin{document}$ u \in C^{m, l}(\bar{ \Omega}) $\end{document}</tex-math></inline-formula> for <inline-formula><tex-math id="M8">\begin{document}$ m \in {\mathbb N} $\end{document}</tex-math></inline-formula> and <inline-formula><tex-math id="M9">\begin{document}$ 0 < l < 1 $\end{document}</tex-math></inline-formula>, our method has an accuracy of <inline-formula><tex-math id="M10">\begin{document}$ {\mathcal O}(h^{\min\{m+l, \, 2\}}) $\end{document}</tex-math></inline-formula> for constant and linear basis functions, while <inline-formula><tex-math id="M11">\begin{document}$ {\mathcal O}(h^{\min\{m+l, \, 4\}}) $\end{document}</tex-math></inline-formula> for quadratic basis functions. Additionally, our method can be readily applied to approximate the generalized fractional Laplacians with symmetric kernel function, and numerical study on the tempered fractional Poisson problem demonstrates its efficiency.</p>