function M3 = intersect_nfsa(M1,M2)
% M3 is the nfsa-intersection of M1 and M2
% see str2nfsa.m for basic nfsa definitions


M3.ALF = unique([M1.ALF M2.ALF]);
nQ1 = length(M1.Q);
nQ2 = length(M2.Q);
M3.Q = sort(coord2ind(list_all_ind_g({M1.Q,M2.Q}),[nQ1 nQ2]))';
M3.q0 = coord2ind([M1.q0 M2.q0],[nQ1 nQ2]);

M3.F = [];
M3 = init_ndelta(M3);
for q1 = 1:nQ1
  for q2 = 1:nQ2
    if (ismember(q1,M1.F)&ismember(q2,M2.F))  % the only difference with union_dfsa
      M3.F(end+1) = coord2ind([q1 q2],[nQ1 nQ2]);
    end
    
    q3 = coord2ind([q1 q2],[nQ1 nQ2]);
    for s = M3.ALF
      R1 = nfsa_delta(M1,q1,s);
      R2 = nfsa_delta(M2,q2,s);
      if (s=='@')
	% every state implicitly goes to itself on epsilon
	R1 = union(R1,q1);
	R2 = union(R2,q2);
      end
      R3 = [];
      for r1=R1
	for r2=R2
	  r3 = coord2ind([r1 r2],[nQ1 nQ2]);
	  R3(end+1) = r3;
	end
      end
      M3 = set_ndelta(M3,q3,s,{R3});
    end
  end
end

M3 = remove_unreachable(M3);
%M3 = minimize_dfsa(M3);

return


function M3 = intersect_nfsa_cheap(M1,M2)
% M3 is the nfsa-intersection of M1 and M2
% see str2nfsa.m for basic nfsa definitions

% simple, obvious, self-explanatory:
% which would be great except negation requires determinization
%notM1 = negate_fsa(M1);
%notM2 = negate_fsa(M2);
%notM3 = union_nfsa(notM1,notM2);
%M3 = negate_fsa(notM3);

d1 = n2dfsa(M1);
d2 = n2dfsa(M2);
d3 = intersect_dfsa(d1,d2);
M3 = d2nfsa(d3);

return