This article presents a new approach to long time wave packet propagation. The methodology relies on energy domain calculations and an on-the-surface straightforward energy to time transformation to provide wave packet time evolution. The adaptive bisection fast Fourier transform method employs selective bisection to create a multiresolution energy grid, dense near resonances. To implement fast Fourier transforms on the nonuniform grid, the uniform grid corresponding to the finest resolution is reconstructed using an iterative interpolation process. By proper choice of the energy grid points, we are able to produce results equivalent to grids of the finest resolution, with far fewer grid points. We have seen savings 20-fold in the number of eigenfunction calculations. Since the method requires computation of energy eigenfunctions, it is best suited for situations where many wave packet propagations are of interest at a fixed small set of points--as in time dependent flux computations. The fast Fourier transform (FFT) algorithm used is an adaptation of the Danielson-Lanczos FFT algorithm to sparse input data. A specific advantage of the adaptive bisection FFT is the possibility of long time wave packet propagations showing slow resonant decay. A method is discussed for obtaining resonance parameters by least squares fitting of energy domain data. The key innovation presented is the means of separating out the smooth background from the sharp resonance structure.