  In this paper, we propose a new class of operator factorization methods to
discretize the integral fractional Laplacian $(-\Delta)^\frac{\alpha}{2}$ for
$\alpha \in (0, 2)$. The main advantage of our method is to easily increase
numerical accuracy by using high-degree Lagrange basis functions, but remain
the scheme structure and computer implementation unchanged. Moreover, our
discretization of the fractional Laplacian results in a symmetric (multilevel)
Toeplitz differentiation matrix, which not only saves memory cost in
simulations but enables efficient computations via the fast Fourier transforms.
The performance of our method in both approximating the fractional Laplacian
and solving the fractional Poisson problems was detailedly examined. It shows
that our method has an optimal accuracy of ${\mathcal O}(h^2)$ for constant or
linear basis functions, while ${\mathcal O}(h^4)$ if quadratic basis functions
are used, with $h$ a small mesh size. Note that this accuracy holds for any
$\alpha \in (0, 2)$ and can be further increased if higher-degree basis
functions are used. If the solution of fractional Poisson problem satisfies $u
\in C^{m, l}(\bar{\Omega})$ for $m \in {\mathbb N}$ and $0 < l < 1$, then our
method has an accuracy of ${\mathcal O}\big(h^{\min\{m+l,\, 2\}}\big)$ for
constant and linear basis functions, while ${\mathcal O}\big(h^{\min\{m+l,\,
4\}}\big)$ for quadratic basis functions. Additionally, our method can be
readily applied to study generalized fractional Laplacians with a symmetric
kernel function, and numerical study on the tempered fractional Poisson problem
demonstrates its efficiency.
