We develop an algorithm for simulating "perfect" random samples from the invariant measure of a Harris recurrent Markov chain. The method uses backward coupling of embedded regeneration times, and works most effectively for finite chains and for stochastically monotone chains even on continuous spaces, where paths may be sandwiched below "upper" and "lower" processes. Examples show that more naive approaches to constructing such bounding processes may be considerably biased, but that the algorithm can be simplified in certain cases to make it easier to run. We give explicit analytic bounds on the backward coupling times in the stochastically monotone case. An application of the simpler algorithm to storage models is given.