We present a fast solver for the Helmholtz equation
Delta u +/- lambda(2) u = f;
in a 3D rectangular box. The method is based on the application of the disc
rete Fourier transform accompanied by a subtraction technique which allows
us to reduce the errors associated with the Gibbs phenomenon and achieve an
y prescribed rate of convergence. The algorithm requires O(N-3 log N) opera
tions, where N is the number of grid points in each direction. We solve a D
irichlet boundary problem for the Helmholtz equation. We also extend the me
thod to the solution of mixed problems, where Dirichlet boundary conditions
are specified on some faces and Neumann boundary conditions are specified
on other faces. High-order accuracy is achieved by a comparatively small nu
mber of points. For example, for the accuracy of 10(-8) the resolution of o
nly 16-32 points in each direction is necessary.