Accurate and efficient numerical simulation of highly nonlinear ultrasound propagation is essential for a wide range of therapeutic and physical ultrasound applications. However, due to large domain sizes and the generation of higher harmonics, such simulations are computationally challenging, particularly in 3-D problems with shock waves. Current numerical methods are based on computationally inefficient uniform meshes that resolve the highest harmonics across the entire spatial domain. To address this challenge, we present an adaptive numerical algorithm for computationally efficient nonlinear acoustic holography. At each propagation step, the algorithm monitors the harmonic content of the acoustic signal and adjusts its discretization parameters accordingly. This enables efficient local resolution of higher harmonics in areas of high nonlinearity while avoiding unnecessary resolution elsewhere. Furthermore, the algorithm actively adapts to the signal's nonlinearity level, eliminating the need for prior reference simulations or information about the spatial distribution of the harmonic content of the acoustic field. The proposed algorithm incorporates an upsampling process in the frequency domain to accommodate the generation of higher harmonics in forward propagation and a downsampling process when higher harmonics are decimated in backward propagation. The efficiency of the algorithm was evaluated for highly nonlinear 3-D problems, demonstrating a significant reduction in computational cost with a nearly 50-fold speedup over a uniform mesh implementation. Our findings enable a more rapid and efficient approach to modeling nonlinear high-intensity focused ultrasound (HIFU) wave propagation.