function L = getcorL(M)
%  L = max {length(simple path) + N_a(n_a+1) + N_b(n_b+1)}
nab = nALF(M);
S = simple_paths(M);
LS = [];

for i=1:length(S)
  s = S{i};
  L = length(s);
  Nab = nloop(M,s);
  L = L + sum(Nab.*(nab+1));
  %L = length(S{i}) + sum(nloop(M,S{i}).*(nab+1));
  LS(end+1) = L;
end
L = max(LS);
return

function Nab = nloop(M,s)
Nab = count_loopsinpath(M,s);
return


function Nab = nloop_old(M)
S = simple_paths(M);
Nab = zeros(size(M.ALF));
for i=1:length(S)
  Nab = max(Nab,count_loopsinpath(M,S{i}));
end
return

function Nab = count_loopsinpath(M,s)
Nab = zeros(size(M.ALF));
q = M.q0;
for i=1:length(s)-1
  s1 = s(i);
  q = dfsa_delta(M,q,s1);
  A = selfloops(M,q);
  ii = ismember(M.ALF,A);
  Nab = Nab + ii;
end
return


function nab = nALF(M)
S = simple_paths(M);
nab = zeros(size(M.ALF));
for i=1:length(S)
  s = S{i};
  [ish,a] = ishomogeneous(s);
  if ish
    j = a-M.ALF(1)+1;
    nab(j) = max(nab(j),length(s));
  end
end

return

function [b,a] = ishomogeneous(s)
if isempty(s)
  a = '';
  b = 0;
else
  a = s(1);
  b = all(s==a);
end
return