  We present a simple discretization scheme for the hypersingular integral
representation of the fractional Laplace operator and solver for the
corresponding fractional Laplacian problem. Through singularity subtraction, we
obtain a regularized integrand that is amenable to the trapezoidal rule with
equispaced nodes, assuming a high degree of regularity in the underlying
function (i.e., $u\in C^6(R^d)$). The resulting quadrature scheme gives a
discrete operator on a regular grid that is translation-invariant and thus can
be applied quickly with the fast Fourier transform. For discretizations of
problems related to space-fractional diffusion on bounded domains, we observe
that the underlying linear system can be efficiently solved via preconditioned
Krylov methods with a preconditioner based on the finite-difference
(non-fractional) Laplacian. We show numerical results illustrating the error of
our simple scheme as well the efficiency of our preconditioning approach, both
for the elliptic (steady-state) fractional diffusion problem and the
time-dependent problem.
