function P = genrandstat(P,m,n)
% generate a random stationary distribution

P = rand(m^n,1); P = P/sum(P);
zz = P*0 + 1e-12; oo = P*0+1;

%f = randn(m^n,1);
f = rand(m^n,1);

Aeq = is_stat(P,m,n);
Aeq = [Aeq;oo'];
Beq = zeros(size(Aeq,1),1); Beq(end)=1;

P = linprog(f,[],[],Aeq,Beq,zz,oo);

return