# on charge d'abord le package noncom,
# qui contient des procedures concernant
# la manipulation des alphabets sous-jacents aux calculs.

with(noncom):

#####################################################
# Given a letter, computes its successor in the 
# underlying alphabet
#####################################################

lyndon[foll_letter] := proc(x)
local i;
   member(x,_ALPHABET,i);
   if i = 'i' then x
   elif i < nops(_ALPHABET) then _ALPHABET[i+1]
   else x
   fi;
end:

#####################################################
# Given a word w and an integer,
# computes its sesquipower, that is the least power
# of the word not exceeding length n, and then
# truncate it to length n.
#####################################################

lyndon[sesqui] := proc(mot,n)
local i, j, mot1;
   if n < length(mot) then mot1 := [substring(mot, 1..n)]
   else
      mot1 := noncom[uncat](mot);
      j := length(mot);
      for i from 1 to n-j do mot1 := [op(mot1), mot1[i]] od;
   fi;
   cat(op(mot1));
end:

#####################################################
# Given a Lyndon word w, computes the Lyndon words
# immediately following it, amongst Lyndon words
# of length <= n.
#####################################################

lyndon[foll_word] := proc(mot, n)
local i, mot1, der;
   mot1 := sesqui(mot, n);
   der := _ALPHABET[nops(_ALPHABET)];
   i := n;
   while substring(mot1,i..i) = der do
      i := i-1;
   od;
   i;
   ``.(substring(mot1,1..i-1)).(lyndon[foll_letter](substring(mot1,i..i)));
end:


#####################################################
# Given an integer n, computes the list of all
# Lyndon words of length at most n, using the
# underlying alphabet _ALPHABET. It follows a
# beautiful work by J.P. Duval, TCS, 1990.
#####################################################

lyndon[generate] := proc(n)
local liste, der, nouveaumot;

   liste := _ALPHABET[1]; der := _ALPHABET[nops(_ALPHABET)];
   nouveaumot := lyndon[foll_word](_ALPHABET[1],n);

while (nouveaumot <> der) do
     liste := liste, nouveaumot;
     nouveaumot := lyndon[foll_word](nouveaumot,n);
     od;
   liste, nouveaumot;
end:

#####################################################
# Factorizes a word into a non-increasing product of
# Lyndon words following J.P. Duval algorithm.
#####################################################

lyndon[factorize] := proc()
local F, w, i, j, k, n;

   w := noncom[uncat](args[1]);
   k := 0;
   n := length(args[1]);
   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
         elif lexorder(w[j], w[i]) then break;
         fi;
      od;
      while k < i do
          F := F, k + (j-i);
          k := k + (j-i);
      od;
   od;
   if nargs = 2 and args[2] = 'N' then
      [F];
   else
      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]);
  
      zip(proc(x, y) cat(op(op(1, x)..op(2, x), y)) end,
          ", [w $'i' = 1..nops([F])]);
   fi   
end:

# lyndon[rstand] := proc()
# local w, s, t;
# 
#    if nargs = 2 then
#       if args[2] = 0 then args[1]
#       else
#          if length(args[1]) = 1 then args[1]
#          else s := lyndon[factorize](substring(args[1], 2..length(args[1])));
#               t := substring(args[1], 1..1);
#               for w in s do
#                  t := [t, lyndon[rstand](w, args[2]-1)];
#               od;
#               t;
#          fi;
#       fi;
#    else
#       if length(args[1]) = 1 then args[1]
#       else s := lyndon[factorize](substring(args[1], 2..length(args[1])));
#            t := substring(args[1], 1..1);
#            for w in s do
#               t := [t, lyndon[rstand](w)];
#            od;
#            t;
#       fi;
#    fi
# end:

lyndon[rstand] := proc(mot)

local i, s, t;

   if length(mot) = 1 then mot
   else s := lyndon[factorize](substring(mot, 2..length(mot)));
        t := substring(mot, 1..1);
        for i from 1 to nops(s) do
           t := [t, lyndon[rstand](op(i, s))];
        od;
        t;
   fi;
end:

lyndon[lstand] := proc(mot)

local i, s, t;

   if length(mot) = 1 then mot
   else s := lyndon[factorize](substring(mot, 1..length(mot)-1));
        t := substring(mot, length(mot)..length(mot));
        for i from nops(s) by -1 to 1 do
           t := [lyndon[lstand](op(i, s)), t];
        od;
        t;
   fi;
end:


######## Definition of help functions associated ########
########       with procedures of the package    ########

`help/text/lyndon` := TEXT(
`HELP FOR: Introduction to the lyndon package`,
``,
`CALLING SEQUENCE:`,
`   <function>(args)`,
`   lyndon[<function>](args)`,
``,
`SYNOPSIS:   `,
`- To use a lyndon function,`,
`  either define that function alone using the command`,
`  with(noncom, <function>), or define all noncom functions`,
`  using the command with(noncom).  Alternatively, invoke`,
`  the function using the long form noncom[<function>].`,
`  This long form notation is necessary whenever there is a`,
`  conflict between a package function name and another function`,
`  used in the same session.`,
``,   
`- The functions available are:`,
``,
`      foll_letter,       sesqui,`,
`      foll_word,            generate`,
``,
`- The package relies on some functionnalities of the noncom`,
`  package for dealing with noncommutative words.`,
`  Mainly, it uses a global variable noncom['_ALPHABET'].`,
`  See noncom['_ALPHABET'] and noncom[alphabet] for more information.`,
``,
`- For more information on a particular function see lyndon[<function>].`,
``,
`- The package deals with Lyndon words. The only relevant function`,
`  to the user is the function generate used to generate all Lyndon`,
`  words up to a given length. A good introduction to the subject is:`,
``,
`     LOTHAIRE M., Combinatorics on Words, Addison-Wesley, 1982.`,
``,
`  More specialized references containing all the ideas in the`,
`  algorithms are:`,
``,
`     DUVAL J.P., Factorizing words over a finite alphabet,`,
`     J. of Algorithms, 1984.`,
``,
`     DUVAL J.P., Generation d'une section des classes de`,
`     conjugaison primitives, TCS, 1990.`,
``
):

`help/lyndon/text/foll_letter` := TEXT(
`HELP FOR: foll_letter, a procedure to obtain the successor,`,
`                           of a letter in the underlying alphabet.`,
``,
`CALLING SEQUENCE:`,
`   foll_letter(x)`,
``,
`PARAMETERS`,
`- x, any letter in the underlying alphabet noncom[_ALPHABET].`,
``,
`SYNOPSIS:`,
`- The function returns the member of the list _ALPHABET`,
`  following the letter x, if it exists. Otherwise, either because`,
`  x does not belong to _ALPHABET or because x is actually the last`,
`  element of _ALPHABET, it returns x itself.`,
``,
`EXAMPLES:`,
`> _ALPHABET;`,
``,
`                                     [a, b]`,
``,
`> foll_letter(a);`,
``,
`                                       b`,
``,
`> foll_letter(b);`,
``,
`                                       b`,
``,
`> foll_letter(c);`,
``,
`                                       c`,
``,
`SEE ALSO: noncom[alphabet], noncom[_ALPHABET].`,
``
):

`help/lyndon/text/sesqui` := TEXT(
`HELP FOR: sesqui, a procedure to compute the sesquipower,`,
`                  of a word up to a given length.`,
``,
`CALLING SEQUENCE:`,
`   sesqui(w, n)`,
``,
`PARAMETERS`,
`- w, any letter in the underlying alphabet noncom[_ALPHABET],`,
`- n, any positive integer.`,
``,
`SYNOPSIS:`,
`- The function computes the n th power of the word w`,
`  and truncate it to length n. It actually check first`,
`  if n is less than the length of w,in which it only`,
`  computes the propoer substring of w.`,
``,
`  Note that this function is independant of the value`,
`  of the underlying alphabet _ALPHABET.`,
``,
`EXAMPLES:`,
`> sesqui(abbab, 2);`,
``,
`                                       ab`,
``,
`> sesqui(abb, 5);  `,
``,
`                                     abbab`,
`> _ALPHABET;`,
``,
``
):

`help/lyndon/text/foll_word` := TEXT(
`HELP FOR: foll_word, a procedure to compute the Lyndon word`,
`                       immediately following a given Lyndon word w`,
`                       (with respect to lexicographical order and`,
`                       amongst words of bounded length n.`,
``,
`CALLING SEQUENCE:`,
`   foll_word(w, n)`,
``,
`PARAMETERS`,
`- w, any noncommutative word.`,
``,
`SYNOPSIS:`,
`- The function computes the Lyndon word immediately following a given`,
`  Lyndon word w (with respect to lexicographical order and amongst words`,
`  of bounded length n. It uses a result of Duval: it suffices to replace,`,
`  in the appropriate sesquipower of w, the rightmost letter not equal to`,
`  the last letter of the alphabet by its immediate successor in the alphabet.`,
``,
`  Note that this function STRONGLY depends on the value of the underlying`,
`  alphabet _ALPHABET.`,
``,
`EXAMPLES:`,
`> foll_word(abb, 5);`,
``,
`                                      abbb`,
``,
`> foll_word(aabb, 5);`,
``,
`                                     aabbb`,
``,
`> foll_word(abbab, 2);`,
``,
`                                       b`,
``
):

`help/lyndon/text/generate` := TEXT(
`HELP FOR: generate, a procedure to compute the list of all Lyndon words`,
`                    of bounded length over the underlying alphabet.`,
``,
`CALLING SEQUENCE:`,
`   generate(n)`,
``,
`PARAMETERS`,
`- n, a positive integer`,
``,
`SYNOPSIS:`,
`- The function computes the list of all Lyndon words of length at most n`,
`  over the underlying alphabet _ALPHABET.`,
`  The underlying algorithms is of J.P. Duval, J. of Algorithms, 1984.`,
``,
`  Note that this function STRONGLY depends on the value of the underlying`,
`  alphabet _ALPHABET.`,
``,
`EXAMPLES:`,
`> _ALPHABET;`,
``,
`                                     [a, b]`,
``,
`> generate(4);`,
``,
`                      a, aaab, aab, aabb, ab, abb, abbb, b`,
``,
`> alphabet([x, z, y]);`,
``,
`                                   [x, y, z]`,
``,
`> generate(3);`,
``,
`          x, xxy, xxz, xy, xyy, xyz, xz, xzy, xzz, y, yyz, yz, yzz, z`,
``
):

save ``.`lyndon.m`;
quit;

C := proc(B)

   map(proc(x) if nops(x) = 1 then lie(x) 
               else noncom[`&t`](lie(op(1, x)), C([op(2..nops(x), x)]))
               fi end, B);
end;