# Get Lyndon factorization of binary vectors of lengths 1, 2, 3, 4, ...
# ans (A211097) will contain the number of factors in the Lyndon factorization 
# ansLs (A211098) will contain the size of the largest (i.e. left-most) factor in the Lyndon factorization
# ansLw (A211099) will contain the largest (i.e. left-most) factor in the Lyndon factorization
# ansRs (A211095 with different offset) will contain the size of the smallest (i.e. right-most) factor in the Lyndon factorization
# ansRw (A211096 with different offset) will contain the smallest (i.e. right-most) factor in the Lyndon factorization

ans:=[];
ansLs:=[];
ansLw:=[];
ansRs:=[];
ansRw:=[];

M:=6;
# study binary vectors of length M
for n from 1 to M do
t1:=binvecs(2^n,2^(n+1)-1);
   for n1 from 1 to nops(t1) do
t2:=t1[n1];
s2:=subsop(1=NULL, t2);
t3:=convert2string(s2);
t4:=factorize(t3);
lprint(n, t2, t4, nops(t4));
ans:=[op(ans),nops(t4)];
ansLs:=[op(ansLs), length(t4[1])];
ansLw:=[op(ansLw), t4[1]];
ansRs:=[op(ansRs), length(t4[nops(t4)])];
ansRw:=[op(ansRw), t4[nops(t4)]];
                             od:
od:

# This uses the following procedures, which of course must be loaded first:


# binvecs Produce binary vectors of all numbers from n1 to n2 (Maple)
binvecs:=proc(n1,n2) local lis,t1,t2,len,n;
lis:=[];
for n from n1 to n2 do
t1:=convert(n,base,2);
len:=nops(t1);
t2:=[seq(t1[len+1-i],i=1..len)];
lis:=[op(lis),t2];
                   od;
lis;
end;


# converts list to string (Maple)
convert2string := proc(x) local a,t1,i;
a:="";
for i from 1 to nops(x) do a:=cat(a,x[i]); od; a; end;
 

# Melancon's "uncat"
uncat := proc (word)
    if length (word) <= 1 then [word]
    else 
         map (proc(x, y ) 
           substring(y, x..x) end,
     [$ 1..length(word)], word)
fi; end:


# Melancon's implementation of Duval's algorithm for Lyndon factorization
factorize := proc(word)
local F, i, j, k, n, w; global uncat;
w := uncat(word); n := length(word); k := 0; F := NULL;

while k  < n do 
   i := k+1;
   j := k+2;
while j < n+1 do
   if lexorder(w[i], w[j]) then
     if w[i] = w [j] then
     i := i+1;
     j := j+1;
else
      i := k+1;
      j := j+1; 
fi 
else break 
fi;
            od;
while k < i do
   F := F, k + (j-i);
   k := k + (j-i);
            od; 
            od;
 
map(proc(x,y) if x < nops (y) then x+1 else nops(y) fi end,
[op(1..nops([F])-1, [F])], w);

zip(proc(x,y) [x, y] end, [1, op(%)], [F]);

map (proc(x,y) cat(op(op(1, x)..op(2,x), y)) end, %, w);
end: