A new algorithm for transferring radiation in three-dimensional space is developed. The algorithm computes radiation sources and sinks using the fast Fourier transform method, based on a formulation in which the integral of any quantity ( such as emissivity or opacity) over any volume may be written in the classic convolution form. The algorithm is fast, with the computational time scaling as N (log N)(2) (where N is the number of grid points of a simulation box), independent of the number of radiation sources. Furthermore, in this formulation one can naturally account for both local radiation sources and diffuse background as well as any extra external sources, in a self-consistent fashion. Finally, the algorithm is completely stable and robust. While the algorithm is generally applicable, we test it thoroughly on an encompassing set of problems relevant for cosmological applications. Tests show that computed results with the present algorithm are, in all cases, in very good agreement with analytic expectations, where they exist. For a realistic cosmological density field the algorithm is in excellent agreement with a high-resolution ray-tracing method. Specifically, radiation flux is guaranteed to propagate in the right direction, and its distribution is shown to be computed accurately. The ionization fronts travel at the correct speeds, with errors no larger than about one cell for all the cases tested. The total number of photons is conserved in the worst case at a similar to10% level and typically at a 1%-5% level over hundreds of time steps. As an added advantage, the accuracy of the results depends weakly on the size of the time step, with a typical cosmological hydrodynamic time step being sufficient.