stationary <- function(trans, eps = 0.005)
{
        n <- nrow(trans)
        if(ncol(trans)!=n)
                stop("transition matrix must be square")
        for(i in 1:n) {
                if(abs(sum(trans[i,  ]) - 1) > eps)
                        stop("each row must sum to one")
        }
# The plan:  need to find a vector pi that statisfies pi %* trans = pi.
# The system we would want to solve is (t(trans) - I) %* pi = 0, where
# pi is now a column vector.  This system is singular, so we will change
# the last row of (t(trans) - I) into a row of one's, and make the right
# hand side vector of zero's into (0,...,0,1).  This last row represents
# the constraint that the pi's sum to one.  We then solve for pi.
# Transpose the transition matrix.
        a <- t(trans)   # subtract one from the diagonal elements
        for(i in 1:n) {
                a[i, i] <- a[i, i] - 1
        }
# change the last row to 1
        a[n,  ] <- 1    # create the right hand side vector
        b <- rep(0, n)
        b[n] <- 1       # solve the system of equations
        x <- solve(a, b)
	x
}