function R = reachable_states(M)
% R contains all those states q s.t. q0 ->* q
ft = fsa_type(M);
switch ft
case 'd'
  R = reachable_states_d(M);
case 'n'
  R = reachable_states_n(M);
end
return


function R = reachable_states_d(M)
% the deterministic version
nQ = length(M.Q);
% we only need to consider paths of length nQ-1
t = 1; b = 1;
Qt = [M.q0];
R = Qt;
for t=1:nQ-1
  Qt1 = [];
  %for q=setdiff(Qt,uR)
  for q=Qt
    %for s=double(M.ALF)
    for s=M.ALF
      r = dfsa_delta(M,q,s);
      Qt1(end+1) = r;
    end
  end
  Qt = unique(Qt1);
  R = unique([R,Qt]);
end
return



function R = reachable_states_n(M)
% the non-deterministic version
nQ = length(M.Q);
% we only need to consider paths of length nQ-1
t = 1; b = 1;
Qt = [M.q0];
R = Qt;
for t=1:nQ-1
  Qt1 = [];
  for q=Qt
    for s=M.ALF
      P = nfsa_delta(M,q,s);
      Qt1 = union(Qt1,P);
    end
  end
  Qt = Qt1;
  R = union(R,Qt);
end
return