We present the theory, algorithm, and numerical experiments for the implementation of a new N-body algorithm, called APM. It scales linearly with the number of particles for the computational effort per time step. The resolution is fully adaptive, with a typical smoothing length comparable to the local interparticle separation. This is accomplished through the use of a dynamical coordinate system, which adjusts itself to the local density distribution. For the Poisson solver a multigrid iteration scheme is used. The scheme is fully data parallel and data local.A specific implementation of this algorithm is described which is available to the astrophysics community. It is optimized for scalar, vector, and shared memory parallel machines. Performance characteristics are compared to traditional particle-mesh (PM) algorithms. The algorithm inherently has a threefold higher computing cost than the PM method with the same number of grid cells, but offers much higher resolution, and currently appears to be amongst the fastest adaptive high-resolution N-body algorithms.