%ever-reaching probability tensor F of a higher order chain
%P is m-th order, n-dimensional transition tensor
%tol is tolerance and is optional, default tol=1e-6
function F=erp(P,tol)
s=size(P);
n=s(1);
m=length(s);
N=n^m;
if nargin<2
    tol=1e-6;
end
F=P;
t=P;
kmax=1000;
for k=1:kmax
    t=bprod(t-diagt(t),P);
    F=F+t;
    if norm(reshape(t,1,N),inf)<tol
        fprintf('F is approximated by the first %3d terms\n',k+1)
        break
    end
end
        
    