The phase shift method (Gazdag, 1978) is based on the solution, in the frequency domain, of an approximation (Claerbout, 1976) to the one‐way wave equation with initial conditions defined by a zero‐offset seismic section. Wave velocity is assumed to be constant within each layer of the section grid and is allowed to vary from layer to layer. Under these conditions, the equation written in the frequency domain reduces to a system of independent ordinary differential equations with initial values that can be solved analytically for each layer. The integration process simply amounts to multiplying the initial values by a complex number of unit modulus. The main advantages of this method are simplicity, stability, and high accuracy, since the precision of the Fourier approximation is limited only by the granularity of the seismic section. From a practical viewpoint of computer implementation, the phase shift method offers a great deal of flexibility. Some accuracy can be traded for speed, as needed, by excluding the waves traveling at an angle that exceeds a specified limit, or by ignoring frequencies at the higher end of the spectrum. The migration of a group of reflections can be excluded as well, since the phase shift angles can be accumulated for application to the frequency representation of the section only within a specified depth interval. The most interesting feature of the method may well reside in the fact that the Fourier coefficients of a seismic section can be treated independently of one another. This is of particular importance if the method is implemented for an array processor. In this case, cross‐sections of the frequency data can be sent to the array processor one‐by‐one (or in groups, as may be suitable) for phase shift operations, at minimal storage cost. A formulation of the method for an array processor is outlined here, which takes advantage of the features mentioned above. It is described in a step‐wise algorithmic fashion, using natural language and standard mathematical notations. Functions for discrete Fourier transformations, evaluation of square roots and circular functions (or fast table look‐up or interpolation), and vector reversal are assumed to exist in the array processor, in addition to elementary vector operations. This formulation can also be used for efficient implementation on conventional computers.