  We develop a novel Monte Carlo parton branching algorithm based on the
Gribov-Levin-Ryskin (GLR) equation. The formulations of both forward evolution
and backward evolution for the GLR equation are presented. The results from the
Monte Carlo implementation of the GLR equation are in full agreement with its
numerical solutions. Our work thus paves the way for developing an event
generator that embodies the saturation effect.
