The overall elasto-plastic behavior of single crystals is governed by individual slips on crystallographic planes, which occur when the resolved shear stress on a critical slip system reaches a certain maximum value. The challenge lies in identifying the activated slip systems for a given load increment since the process involves selection from a pool of linearly dependent slip systems. In this paper, we use an “ultimate algorithm” for the numerical integration of the elasto-plastic constitutive equation for single crystals. The term ultimate indicates exact integration of the elasto-plastic constitutive equation and explicit tracking of the sequence of slip system activation. We implement the algorithm into a finite element code and report the performance for polycrystals subjected to complicated loading paths including non-proportional and reverse/cyclic loading at different crystal orientations. It is shown that the ultimate algorithm is comparable to the widely used radial return algorithm for J2 plasticity in terms of global numerical stability.