In this paper, we propose a new class of operator factorization methods to discretize the integral fractional Laplacian (−∆)α/2 for α ∈ (0, 2). 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 O(h2), while O(h4) for quadratic basis functions with h a small mesh size. This accuracy can be achieved for any α ∈ (0, 2) 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 u ∈ Cm,l((equation presented)) for m ∈ ℕ and 0 < l < 1, our method has an accuracy of O(hmin{m+l, 2}) for constant and linear basis functions, while O(hmin{m+l, 4}) 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.