# Get Lyndon factorization of binary expansion of numbers from 0 to M
# ans (A211100) will contain the number of factors in the Lyndon factorization of n
# ansLs (always 1) will contain the size of the largest (i.e. left-most) factor in the Lyndon factorization of n
# ansRs (A211095) will contain the size of the smallest (i.e. right-most) factor in the Lyndon factorization of n
# ansRw (A211096) will contain the smallest (i.e. right-most) factor in the Lyndon factorization of n

ans:=[];
ansLs:=[];
ansRs:=[];
ansRw:=[];
M:=200;
for n from 0 to M do
t1:=binvecs(n,n);
t2:=op(t1);
t3:=convert2string(t2);
t4:=factorize(t3);
#lprint(n, t1, t4, nops(t4));
ans:=[op(ans),nops(t4)];
ansLs:=[op(ansLs), length(t4[1])];
ansRs:=[op(ansRs), length(t4[nops(t4)])];
ansRw:=[op(ansRw), t4[nops(t4)]];
od:

# This uses the following procedures


# 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: