  We introduce a general purpose algorithm for rapidly computing certain types
of oscillatory integrals which frequently arise in problems connected to wave
propagation and general hyperbolic equations. The problem is to evaluate
numerically a so-called Fourier integral operator (FIO) of the form $\int
e^{2\pi i \Phi(x,\xi)} a(x,\xi) \hat{f}(\xi) \mathrm{d}\xi$ at points given on
a Cartesian grid. Here, $\xi$ is a frequency variable, $\hat f(\xi)$ is the
Fourier transform of the input $f$, $a(x,\xi)$ is an amplitude and
$\Phi(x,\xi)$ is a phase function, which is typically as large as $|\xi|$;
hence the integral is highly oscillatory at high frequencies. Because an FIO is
a dense matrix, a naive matrix vector product with an input given on a
Cartesian grid of size $N$ by $N$ would require $O(N^4)$ operations.
  This paper develops a new numerical algorithm which requires $O(N^{2.5} \log
N)$ operations, and as low as $O(\sqrt{N})$ in storage space. It operates by
localizing the integral over polar wedges with small angular aperture in the
frequency plane. On each wedge, the algorithm factorizes the kernel $e^{2 \pi i
\Phi(x,\xi)} a(x,\xi)$ into two components: 1) a diffeomorphism which is
handled by means of a nonuniform FFT and 2) a residual factor which is handled
by numerical separation of the spatial and frequency variables. The key to the
complexity and accuracy estimates is that the separation rank of the residual
kernel is \emph{provably independent of the problem size}. Several numerical
examples demonstrate the efficiency and accuracy of the proposed methodology.
We also discuss the potential of our ideas for various applications such as
reflection seismology.
