%
%   隨機數學(30250143, 07Fall) 大作業程序
%
%   作者： 王  音 2004011448 
%         葛元愷 2004011425 <geyuank@gmail.com>
%
%   All rights reserved. 2008.1

%   參看文獻《馬爾可夫基礎及其應用》P151-P161

function result = Markov_Return(P)

[m n] = size(P);

if m~=n,
    error('P must be a square-matrix (m=n)');
end

N = n;
M = 0;
Z = calc_Z(P);
Z0= Z;

while(2^M < (N-1)*(N-1)),   %計算出符合要求的最小Z^(2^M)
    Z = calc_Z(Z*Z);
    M = M + 1;
end

disp(' ');              %空一行
Diff = 1:n;

while (sum(sum(Z))>0),  %判斷Z是否已經為全零陣

    L=Z;
    
    C = find(sum(Z)==max(sum(Z)));
    C0 = C(1);          %找到一個常返態
    d = 1;              %記錄周期信息
    
    G0 = calc_G(L,C0);  
    C  = G0;

    L = calc_Z(L*Z0);
    Gd = calc_G(L,C0);

    while Gd~=G0,   %計算各個Gd，參看P152中間
        C = union(C,Gd);
        L = calc_Z(L*Z0);
        Gd = calc_G(L,C0);
        d = d+1;
    end

    D = C;          %計算D，參看P153上側
    for i=1:length(C),
        for j=1:n,
            if Z(j,C(i))==1,
                D = union(D,j);
            end
        end
    end

    for i=1:length(D),%把D中對于狀態的行和列都置0，見P153中間
        Z(:,D(i))=0;
        Z(D(i),:)=0;
    end
    
    disp(sprintf('找到常返態閉集，周期為%d',d));    
    disp(C);
        
    Diff = setdiff(Diff,C); %剩下的狀態

end

if Diff,                    %打印出剩下的非常返的狀態
    disp('其余非常返閉集');
    disp(Diff);
end