--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

import OrderedSet
import RetractableTo
)abbrev domain MAGMA Magma
++ Author: Michel Petitot (petitot@lifl.fr).
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This type is the basic representation of 
++ parenthesized words (binary trees over arbitrary symbols)
++ useful in \spadtype{LiePolynomial}. \newline Author: Michel Petitot (petitot@lifl.fr).

Magma(VarSet:OrderedSet):Public == Private where
   WORD ==> OrderedFreeMonoid(VarSet)
   EX   ==> OutputForm

   Public == Join(OrderedSet,RetractableTo VarSet) with
      *           : ($,$) -> $
        ++ \axiom{x*y} returns the tree \axiom{[x,y]}.
      coerce        : $ -> WORD
        ++ \axiom{coerce(x)} returns the element of \axiomType{OrderedFreeMonoid}(VarSet) 
        ++ corresponding to \axiom{x} by removing parentheses.
      first         : $ -> VarSet
        ++ \axiom{first(x)} returns the first entry of the tree \axiom{x}.
      left          : $ -> $
        ++ \axiom{left(x)} returns left subtree of \axiom{x} or
        ++ error if \axiomOpFrom{retractable?}{Magma}(\axiom{x}) is true.
      length        : $ -> PositiveInteger
        ++ \axiom{length(x)} returns the number of entries in \axiom{x}.
      lexico        : ($,$) -> Boolean
        ++ \axiom{lexico(x,y)} returns \axiom{true} iff  \axiom{x} is smaller than 
        ++ \axiom{y} w.r.t. the lexicographical ordering induced by \axiom{VarSet}. 
        ++ N.B. This operation does not take into account the tree structure of
        ++ its arguments. Thus this is not a total ordering.
      mirror        : $ -> $
        ++ \axiom{mirror(x)} returns the reversed word of \axiom{x}. 
        ++ That is \axiom{x} itself if \axiomOpFrom{retractable?}{Magma}(\axiom{x}) is true and
        ++ \axiom{mirror(z) * mirror(y)} if \axiom{x} is \axiom{y*z}.
      rest          : $ -> $
        ++ \axiom{rest(x)} return \axiom{x} without the first entry or 
        ++ error if \axiomOpFrom{retractable?}{Magma}(\axiom{x}) is true.
      retractable?  : $ -> Boolean
        ++ \axiom{retractable?(x)} tests if \axiom{x} is a tree with only one entry.
      right         : $ -> $
        ++ \axiom{right(x)} returns right subtree of \axiom{x} or 
        ++ error if \axiomOpFrom{retractable?}{Magma}(\axiom{x}) is true.
      varList       : $ -> List VarSet
        ++ \axiom{varList(x)} returns the list of distinct entries of \axiom{x}.

   Private == add
    -- representation
      VWORD := Record(left:$ ,right:$)
      Rep:= Union(VarSet,VWORD)  

      recursif: ($,$) -> Boolean

    -- define
      x:$ = y:$ ==
        x case VarSet => 
           y case VarSet => x::VarSet = y::VarSet
           false
        y case VWORD => x::VWORD = y::VWORD
        false
 
      varList x == 
        x case VarSet => [x::VarSet]
        lv: List VarSet := setUnion(varList x.left, varList x.right)
        sort!(lv)

      left x == 
        x case VarSet => error "x has only one entry"
        x.left

      right x == 
        x case VarSet => error "x has only one entry"
        x.right
      retractable? x == (x case VarSet)

      retract x ==
         x case VarSet => x::VarSet
         error "Not retractable"

      retractIfCan x == (retractable? x => x::VarSet ; "failed")
      coerce(l:VarSet):$  == l

      mirror x ==
        x case VarSet => x
        [mirror x.right, mirror x.left]$VWORD

      coerce(x:$): WORD ==
        x case VarSet => x::VarSet::WORD
        x.left::WORD * x.right::WORD

      coerce(x:$):EX ==
         x case VarSet => x::VarSet::EX
         bracket [x.left::EX, x.right::EX]

      x * y == [x,y]$VWORD

      first x ==
         x case VarSet => x::VarSet
         first x.left

      rest x ==
         x case VarSet => error "rest$Magma: inexistant rest"
         lx:$ := x.left
         lx case VarSet => x.right
         [rest lx , x.right]$VWORD

      length x ==
         x case VarSet => 1
         length(x.left) + length(x.right)

      recursif(x,y) ==    
         x case VarSet => 
            y case VarSet => x::VarSet < y::VarSet
            true
         y case VarSet => false
         x.left =  y.left =>  x.right <  y.right
         x.left < y.left

      lexico(x,y) ==      	-- peut etre amelioree !!!!!!!!!!!
         x case VarSet => 
            y case VarSet => x::VarSet < y::VarSet
            x::VarSet <= first y
         y case VarSet => first x < retract y
         fx:VarSet := first x ; fy:VarSet := first y 
         fx = fy => lexico(rest x , rest y)
         fx < fy 

      x < y ==                 	-- recursif par longueur
         lx,ly: PositiveInteger
         lx:= length x ; ly:= length y
         lx = ly => recursif(x,y)
         lx < ly 

import OrderedSet
import RetractableTo
import Boolean
import Magma
)abbrev domain LWORD LyndonWord
++ Author: Michel Petitot (petitot@lifl.fr).
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References: Free Lie Algebras by C. Reutenauer (Oxford science publications).
++ Description:
++ Lyndon words over arbitrary (ordered) symbols:
++ see Free Lie Algebras by C. Reutenauer (Oxford science publications).
++ A Lyndon word is a word which is smaller than any of its right factors
++ w.r.t. the pure lexicographical ordering.
++ If \axiom{a} and \axiom{b} are two Lyndon words such that \axiom{a < b}
++ holds w.r.t lexicographical ordering then \axiom{a*b} is a Lyndon word.
++ Parenthesized Lyndon words can be generated from symbols by using the following
++ rule: \axiom{[[a,b],c]} is a Lyndon word iff \axiom{a*b < c <= b} holds.
++ Lyndon words are internally represented by binary trees using the
++ \spadtype{Magma} domain constructor.
++ Two ordering are provided: lexicographic and 
++ length-lexicographic. \newline 
++ Author : Michel Petitot (petitot@lifl.fr).

LyndonWord(VarSet:OrderedSet):Public == Private where
   OFMON ==> OrderedFreeMonoid(VarSet)
   PI    ==> PositiveInteger
   NNI   ==> NonNegativeInteger
   I     ==> Integer
   OF    ==> OutputForm
   ARRAY1==> OneDimensionalArray

   Public == Join(OrderedSet,RetractableTo VarSet) with
      retractable?  : $ -> Boolean
        ++ \axiom{retractable?(x)} tests if \axiom{x} is a tree with only one entry.
      left          : $ -> $
        ++ \axiom{left(x)} returns left subtree of \axiom{x} or
        ++ error if \axiomOpFrom{retractable?}{LyndonWord}(\axiom{x}) is true.
      right  :  $ -> $
        ++ \axiom{right(x)} returns right subtree of \axiom{x} or
        ++ error if \axiomOpFrom{retractable?}{LyndonWord}(\axiom{x}) is true.
      length :  $ -> PI
        ++ \axiom{length(x)} returns the number of entries in \axiom{x}.
      lexico :  ($,$) -> Boolean 
        ++ \axiom{lexico(x,y)} returns \axiom{true} iff  \axiom{x} is smaller than 
        ++ \axiom{y} w.r.t. the lexicographical ordering induced by \axiom{VarSet}. 
      coerce :  $ -> OFMON
        ++ \axiom{coerce(x)} returns the element of \axiomType{OrderedFreeMonoid}(VarSet) 
        ++ corresponding to \axiom{x}.
      coerce :  $ -> Magma VarSet
        ++ \axiom{coerce(x)} returns the element of \axiomType{Magma}(VarSet)
        ++ corresponding to \axiom{x}.
      factor :  OFMON -> List $  
        ++ \axiom{factor(x)} returns the decreasing factorization into Lyndon words. 
      lyndon?:  OFMON -> Boolean
        ++ \axiom{lyndon?(w)} test if \axiom{w} is a Lyndon word.
      lyndon :  OFMON -> $
        ++ \axiom{lyndon(w)} convert \axiom{w} into a Lyndon word, 
        ++ error if \axiom{w} is not a Lyndon word.
      lyndonIfCan : OFMON -> Union($, "failed")
        ++ \axiom{lyndonIfCan(w)} convert \axiom{w} into a Lyndon word.
      varList     : $ -> List VarSet
        ++ \axiom{varList(x)} returns the list of distinct entries of \axiom{x}.
      LyndonWordsList1: (List VarSet, PI)  -> ARRAY1 List $
        ++ \axiom{LyndonWordsList1(vl, n)} returns an array of lists of Lyndon
        ++ words over the alphabet \axiom{vl}, up to order \axiom{n}.
      LyndonWordsList : (List VarSet, PI)  -> List $
        ++ \axiom{LyndonWordsList(vl, n)} returns the list of Lyndon
        ++ words over the alphabet \axiom{vl}, up to order \axiom{n}.

   Private == Magma(VarSet) add
     -- Representation
       Rep:= Magma(VarSet)

     -- Fonctions locales
       LetterList : OFMON -> List VarSet
       factor1    : (List $, $, List $) -> List $

     -- Definitions
       lyndon? w ==
         w = 1$OFMON => false
         f: OFMON := rest w
         while not one? f repeat
           not lexico(w,f) => return false
           f := rest f
         true

       lyndonIfCan w ==
         l: List $ := factor w
         # l = 1 => first l
         "failed"

       lyndon w ==
         l: List $ := factor w
         # l = 1 => first l
         error "This word is not a Lyndon word"

       LetterList w ==
         w = 1 => []
         cons(first w , LetterList rest w)

       factor1 (gauche, x, droite) == 
         g: List $ := gauche; d: List $ := droite
         while not null g repeat             ++ (l in g or l=x) et u in d 
           lexico(  g.first , x ) =>         ++  => right(l) >= u 
              x  := g.first *$Rep x          -- crochetage
              null(d) => g := rest g
              g := cons( x, rest g)          -- mouvement a droite
              x  := first d
              d := rest d
           d := cons( x , d)                 -- mouvement a gauche
           x  := first g
           g := rest g
         return cons(x, d)

       factor w ==
         w = 1 => []
         l : List $ := reverse [ u::$ for u in LetterList w]
         factor1( rest l, first l , [] )
      
       x < y ==                     -- lexicographique par longueur
         lx,ly: PI
         lx:= length x ; ly:= length y
         lx = ly => lexico(x,y)
         lx < ly
 
       coerce(x:$):OF == bracket(x::OFMON::OF)
       coerce(x:$):Magma VarSet == x::Rep

       LyndonWordsList1 (vl,n) ==    -- a ameliorer !!!!!!!!!!!
            null vl => error "empty list"
            base: ARRAY1 List $ := new(n::I::NNI ,[])
           
           -- mots de longueur 1
            lbase1:List $ := [w::$ for w in sort(vl)]
            base.1 := lbase1

           -- calcul des mots de longueur ll
            for ll in 2..n:I  repeat 
               lbase1 := []   
               for a in base(1) repeat              -- lettre + mot
                  for b in base(ll-1) repeat
                     if lexico(a,b) then lbase1:=cons(a*b,lbase1)

               for i in 2..ll-1 repeat              -- mot + mot
                 for a in base(i) repeat             
                   for b in base(ll-i) repeat
                     if lexico(a,b) and (lexico(b,right a) or b = right a ) 
                     then lbase1:=cons(a*b,lbase1)
 
               base(ll):= sort!(lexico, lbase1)
            return base
           
       LyndonWordsList (vl,n) ==
           v:ARRAY1 List $ := LyndonWordsList1(vl,n)
           "append"/ [v.i for i in 1..n] 

import CommutativeRing
import Field
)abbrev category LIECAT LieAlgebra
++ Author: Michel Petitot (petitot@lifl.fr).
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ The category of Lie Algebras.
++ It is used by the following domains of non-commutative algebra:
++ \axiomType{LiePolynomial} and 
++ \axiomType{XPBWPolynomial}. \newline Author : Michel Petitot (petitot@lifl.fr).
LieAlgebra(R: CommutativeRing): Category ==  Module(R) with
  --attributes
    NullSquare 
          ++ \axiom{NullSquare} means that \axiom{[x,x] = 0} holds.
    JacobiIdentity 
          ++ \axiom{JacobiIdentity} means that \axiom{[x,[y,z]]+[y,[z,x]]+[z,[x,y]] = 0} holds.
  --exports
    construct:  ($,$) -> $
          ++ \axiom{construct(x,y)} returns the Lie bracket of \axiom{x} and \axiom{y}.
    if R has Field then 
       /   :  ($,R) -> $
         ++ \axiom{x/r} returns the division of \axiom{x} by \axiom{r}.


  add
    if R has Field then x / r == inv(r)$R * x

import OrderedSet
import CommutativeRing
import LieAlgebra
)abbrev category FLALG FreeLieAlgebra
++ Author: Michel Petitot (petitot@lifl.fr)
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ The category of free Lie algebras.
++ It is used by domains of non-commutative algebra:
++ \spadtype{LiePolynomial} and 
++ \spadtype{XPBWPolynomial}. \newline Author: Michel Petitot (petitot@lifl.fr)

FreeLieAlgebra(VarSet:OrderedSet, R:CommutativeRing) :Category == CatDef where
   XRPOLY  ==> XRecursivePolynomial(VarSet,R)
   XDPOLY  ==> XDistributedPolynomial(VarSet,R)
   RN      ==> Fraction Integer
   LWORD   ==> LyndonWord(VarSet)

   CatDef ==  Join(LieAlgebra(R)) with
      coef      : (XRPOLY , $) -> R
         ++ \axiom{coef(x,y)} returns the scalar product of \axiom{x} by \axiom{y},
         ++ the set of words being regarded as an orthogonal basis.
      coerce    : VarSet -> $
         ++ \axiom{coerce(x)} returns \axiom{x} as a Lie polynomial.
      coerce    : $ -> XDPOLY
         ++ \axiom{coerce(x)} returns \axiom{x} as distributed polynomial.
      coerce    : $ -> XRPOLY 
         ++ \axiom{coerce(x)} returns \axiom{x} as a recursive polynomial.
      degree    : $ -> NonNegativeInteger
         ++ \axiom{degree(x)} returns the greatest length of a word in the support of \axiom{x}.
      --if R has Module(RN) then
      --  Hausdorff : ($,$,PositiveInteger) -> $
      lquo      : (XRPOLY , $) -> XRPOLY
         ++ \axiom{lquo(x,y)} returns the left simplification of \axiom{x} by \axiom{y}.
      rquo      : (XRPOLY , $) -> XRPOLY
         ++ \axiom{rquo(x,y)} returns the right simplification of \axiom{x} by \axiom{y}.
      LiePoly   : LWORD -> $
         ++ \axiom{LiePoly(l)} returns the bracketed form of \axiom{l} as a Lie polynomial.
      mirror    : $ -> $
         ++ \axiom{mirror(x)} returns \axiom{Sum(r_i mirror(w_i))}
         ++ if \axiom{x} is \axiom{Sum(r_i w_i)}.
      trunc     : ($, NonNegativeInteger) -> $
         ++ \axiom{trunc(p,n)} returns the polynomial \axiom{p} 
         ++ truncated at order \axiom{n}.
      varList   : $ -> List VarSet
         ++ \axiom{varList(x)} returns the list of distinct entries of \axiom{x}.
      eval      : ($, VarSet, $) -> $
         ++ \axiom{eval(p, x, v)} replaces \axiom{x} by \axiom{v}  in \axiom{p}.
      eval      : ($, List VarSet, List $) -> $
         ++ \axiom{eval(p, [x1,...,xn], [v1,...,vn])} replaces \axiom{xi} by \axiom{vi}
         ++ in \axiom{p}.

import OrderedSet
import XPolynomialsCat
import NonNegativeInteger
)abbrev package XEXPPKG XExponentialPackage
++ Author: Michel Petitot (petitot@lifl.fr).
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package provides computations of logarithms and exponentials 
++ for polynomials in non-commutative 
++ variables. \newline Author: Michel Petitot (petitot@lifl.fr).

XExponentialPackage(R, VarSet, XPOLY): Public == Private where
    RN     ==> Fraction Integer 
    NNI    ==> NonNegativeInteger
    I      ==> Integer
    R      : Join(Ring, Module RN)
    -- R      : Field
    VarSet : OrderedSet
    XPOLY  : XPolynomialsCat(VarSet, R)

    Public == with
       exp:  (XPOLY, NNI) -> XPOLY
         ++ \axiom{exp(p, n)} returns the exponential of \axiom{p}
         ++ truncated at order \axiom{n}.
       log:  (XPOLY, NNI) -> XPOLY
         ++ \axiom{log(p, n)} returns the logarithm of \axiom{p}
         ++ truncated at order \axiom{n}.
       Hausdorff: (XPOLY, XPOLY, NNI) -> XPOLY
         ++ \axiom{Hausdorff(a,b,n)} returns log(exp(a)*exp(b))
         ++ truncated at order \axiom{n}.

    Private == add
  
        log (p,n) ==
           p1 : XPOLY := p - 1
           not quasiRegular? p1 => 
             error "constant term <> 1, impossible log"
           s : XPOLY := 0       -- resultat
           k : I := n :: I 
           for i in 1 .. n repeat
              k1 :RN := 1/k
              k2 : R := k1 * 1$R
              s := trunc( trunc(p1,i) * (k2 :: XPOLY - s) , i)
              k := k - 1
           s

        exp (p,n) ==
           not quasiRegular? p => 
             error "constant term <> 0, exp impossible"
           p = 0 => 1
           s : XPOLY := 1$XPOLY       -- resultat
           k : I := n :: I
           for i in 1 .. n repeat
              k1 :RN := 1/k
              k2 : R := k1 * 1$R
              s := trunc( 1 +$XPOLY k2 * trunc(p,i) * s , i)
              k := k - 1
           s

        Hausdorff(p,q,n) ==
           p1: XPOLY := exp(p,n)
           q1: XPOLY := exp(q,n)
           log(p1*q1, n)

import OrderedSet
import CommutativeRing
import FreeLieAlgebra
import FreeModuleCat
import LyndonWord
)abbrev domain LPOLY LiePolynomial
++ Author: Michel Petitot (petitot@lifl.fr).
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:Free Lie Algebras by C. Reutenauer (Oxford science publications). 
++ Description:
++ This type supports Lie polynomials in Lyndon basis
++ see Free Lie Algebras by C. Reutenauer 
++ (Oxford science publications). \newline Author: Michel Petitot (petitot@lifl.fr).

LiePolynomial(VarSet:OrderedSet, R:CommutativeRing) : Public == Private where
   MAGMA   ==> Magma(VarSet)
   LWORD   ==> LyndonWord(VarSet)
   WORD    ==> OrderedFreeMonoid(VarSet)
   XDPOLY  ==> XDistributedPolynomial(VarSet,R)
   XRPOLY  ==> XRecursivePolynomial(VarSet,R)
   NNI     ==> NonNegativeInteger
   RN      ==> Fraction Integer
   EX      ==> OutputForm
   TERM    ==> Record(k: LWORD, c: R)

   Public == Join(FreeLieAlgebra(VarSet,R), FreeModuleCat(R,LWORD)) with
      LiePolyIfCan: XDPOLY -> Union($, "failed")
        ++ \axiom{LiePolyIfCan(p)} returns \axiom{p} in Lyndon basis
        ++ if \axiom{p} is a Lie polynomial, otherwise \axiom{"failed"}
        ++ is returned.
      construct: (LWORD, LWORD) -> $
        ++ \axiom{construct(x,y)} returns the Lie bracket \axiom{[x,y]}.
      construct: (LWORD, $) -> $
        ++ \axiom{construct(x,y)} returns the Lie bracket \axiom{[x,y]}.
      construct: ($, LWORD) -> $     
        ++ \axiom{construct(x,y)} returns the Lie bracket \axiom{[x,y]}.

   Private ==  FreeModule1(R, LWORD) add       
        import(TERM)

      --representation
        Rep :=  List TERM

      -- fonctions locales
        cr1 : (LWORD, $    ) -> $
        cr2 : ($, LWORD    ) -> $
        crw : (LWORD, LWORD) -> $     -- crochet de 2 mots de Lyndon
        DPoly: LWORD -> XDPOLY
        lquo1: (XRPOLY , LWORD) -> XRPOLY
        lyndon: (LWORD, LWORD) -> $
        makeLyndon: (LWORD, LWORD) -> LWORD
        rquo1: (XRPOLY , LWORD) -> XRPOLY
        RPoly: LWORD -> XRPOLY
        eval1: (LWORD, VarSet, $) -> $                     -- 08/03/98
        eval2: (LWORD, List VarSet, List $) -> $           -- 08/03/98


      -- Evaluation
        eval1(lw,v,nv) ==                                  -- 08/03/98
          not member?(v, varList(lw)$LWORD) => LiePoly lw
          (s := retractIfCan(lw)$LWORD) case VarSet => 
             if (s::VarSet) = v then nv else LiePoly lw 
          l: LWORD := left lw
          r: LWORD := right lw
          construct(eval1(l,v,nv), eval1(r,v,nv))

        eval2(lw,lv,lnv) ==                                -- 08/03/98
          p: Integer
          (s := retractIfCan(lw)$LWORD) case VarSet =>
             p := position(s::VarSet, lv)$List(VarSet) 
             if p=0 then lw::$ else elt(lnv,p)$List($)
          l: LWORD := left lw
          r: LWORD := right lw
          construct(eval2(l,lv,lnv), eval2(r,lv,lnv))

        eval(p:$, v: VarSet, nv: $): $ ==                  -- 08/03/98
          +/ [t.c * eval1(t.k, v, nv) for t in p]

        eval(p:$, lv: List(VarSet), lnv: List($)): $ ==    -- 08/03/98
          +/ [t.c * eval2(t.k, lv, lnv) for t in p]

        lquo1(p,lw) ==
          constant? p => 0$XRPOLY
          retractable? lw => lquo(p, retract lw)$XRPOLY
          lquo1(lquo1(p, left lw),right lw) - lquo1(lquo1(p, right lw),left lw)  
        rquo1(p,lw) ==
          constant? p => 0$XRPOLY
          retractable? lw => rquo(p, retract lw)$XRPOLY
          rquo1(rquo1(p, left lw),right lw) - rquo1(rquo1(p, right lw),left lw)

        coef(p, lp) == coef(p, lp::XRPOLY)$XRPOLY

        lquo(p, lp) ==
          lp = 0 => 0$XRPOLY
          +/ [t.c * lquo1(p,t.k) for t in lp]
 
        rquo(p, lp) ==
          lp = 0 => 0$XRPOLY
          +/ [t.c * rquo1(p,t.k) for t in lp] 

        LiePolyIfCan p ==         -- inefficace a cause de la rep. de XDPOLY
           not quasiRegular? p => "failed"
           p1: XDPOLY := p ; r:$ := 0
           while p1 ~= 0 repeat
             t: Record(k:WORD, c:R) := mindegTerm p1
             w: WORD := t.k; coef:R := t.c
             (l := lyndonIfCan(w)$LWORD) case "failed" => return "failed"
             lp:$ := coef * LiePoly(l::LWORD)
             r := r + lp 
             p1 := p1 - lp::XDPOLY 
           r
 
      --definitions locales
        makeLyndon(u,v) == (u::MAGMA * v::MAGMA) pretend LWORD
 
        crw(u,v) ==               -- u et v sont des mots de Lyndon
          u = v => 0
          lexico(u,v) => lyndon(u,v)
          - lyndon (v,u)

        lyndon(u,v) ==            -- u et v sont des mots de Lyndon tq u < v
          retractable? u => monom(makeLyndon(u,v),1)
          u1: LWORD := left u
          u2: LWORD := right u
          lexico(u2,v) => cr1(u1, lyndon(u2,v)) + cr2(lyndon(u1,v), u2)
          monom(makeLyndon(u,v),1)
           
        cr1 (l, p) ==
            +/[t.c * crw(l, t.k) for t in p]

        cr2 (p, l) ==
            +/[t.c * crw(t.k, l) for t in p]

        DPoly w ==
           retractable? w => retract(w) :: XDPOLY 
           l:XDPOLY := DPoly left w
           r:XDPOLY := DPoly right w
           l*r - r*l

        RPoly w ==
           retractable? w => retract(w) :: XRPOLY 
           l:XRPOLY := RPoly left w
           r:XRPOLY := RPoly right w
           l*r - r*l 
    
      -- definitions

        coerce(v:VarSet) == monom(v::LWORD , 1)

        construct(x:$ , y:$):$ ==
            +/[t.c * cr1(t.k, y) for t in x]

        construct(l:LWORD , p:$):$ == cr1(l,p) 
        construct(p:$ , l:LWORD):$ == cr2(p,l)
        construct(u:LWORD , v:LWORD):$ == crw(u,v)

        coerce(p:$):XDPOLY ==
            +/ [t.c * DPoly(t.k) for t in p]

        coerce(p:$):XRPOLY ==
            +/ [t.c * RPoly(t.k) for t in p]

        LiePoly(l) == monom(l,1)

        varList p ==
          le : List VarSet := "setUnion"/[varList(t.k)$LWORD for t in p]
          sort(le)$List(VarSet)

        mirror p ==
          [[t.k, (odd? length t.k => t.c; -t.c)]$TERM for t in p]

        trunc(p, n) ==
          degree(p) > n => trunc( reductum p , n)
          p

        degree p == 
          null p => 0
          length( p.first.k)$LWORD

      --  ListOfTerms p == p pretend List TERM
        
--        coerce(x) : EX ==
--           null x => (0$R) :: EX
--           le : List EX := nil
--           for rec in x repeat
--             rec.c = 1$R => le := cons(rec.k :: EX, le)
--             le := cons(mkBinary("*"::EX,  rec.c :: EX, rec.k :: EX), le)
--           1 = #le => first le
--           mkNary("+" :: EX,le)

import OrderedSet
import RetractableTo
import LyndonWord
)abbrev domain PBWLB PoincareBirkhoffWittLyndonBasis
++ Author: Michel Petitot (petitot@lifl.fr).
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This domain provides the internal representation
++ of polynomials in non-commutative variables written
++ over the Poincare-Birkhoff-Witt basis.
++ See the \spadtype{XPBWPolynomial} domain constructor.
++ See Free Lie Algebras by C. Reutenauer 
++ (Oxford science publications). \newline Author: Michel Petitot (petitot@lifl.fr).

PoincareBirkhoffWittLyndonBasis(VarSet: OrderedSet): Public == Private where
   WORD    ==> OrderedFreeMonoid(VarSet)
   LWORD   ==> LyndonWord(VarSet)
   LWORDS  ==> List(LWORD)
   PI      ==> PositiveInteger
   NNI     ==> NonNegativeInteger
   EX      ==> OutputForm

   Public == Join(OrderedSet, RetractableTo LWORD) with
      1: constant -> %
         ++ \spad{1} returns the empty list.
      coerce       : $ -> WORD
         ++ \spad{coerce([l1]*[l2]*...[ln])} returns the word \spad{l1*l2*...*ln},
         ++ where \spad{[l_i]} is the backeted form of the Lyndon word \spad{l_i}.
      coerce       : VarSet -> $
         ++ \spad{coerce(v)} return \spad{v}
      first        : $ -> LWORD
         ++ \spad{first([l1]*[l2]*...[ln])} returns the Lyndon word \spad{l1}.
      length       : $ -> NNI
         ++ \spad{length([l1]*[l2]*...[ln])} returns the length of the word \spad{l1*l2*...*ln}.
      ListOfTerms  : $ -> LWORDS
         ++ \spad{ListOfTerms([l1]*[l2]*...[ln])} returns the list of words \spad{l1, l2, .... ln}.
      rest         : $ -> $
         ++ \spad{rest([l1]*[l2]*...[ln])} returns the list \spad{l2, .... ln}.
      retractable? : $ -> Boolean
         ++ \spad{retractable?([l1]*[l2]*...[ln])} returns true iff \spad{n}  equals \spad{1}.
      varList      : $ -> List VarSet
         ++ \spad{varList([l1]*[l2]*...[ln])} returns the list of
         ++ variables in the word \spad{l1*l2*...*ln}.
   
   Private == add

    -- Representation
     Rep := LWORDS

    -- Locales
     recursif: ($,$) -> Boolean

    -- Define
     1 == nil

     x = y == x =$Rep y

     varList x ==
        null x => nil
        le: List VarSet := "setUnion"/ [varList$LWORD l for l in x]

     first x == first(x)$Rep
     rest x == rest(x)$Rep

     coerce(v: VarSet):$ == [ v::LWORD ]
     coerce(l: LWORD):$ == [l]
     ListOfTerms(x:$):LWORDS == x pretend LWORDS      

     coerce(x:$):WORD ==
       null x => 1
       x.first :: WORD *$WORD coerce(x.rest)

     coerce(x:$):EX ==
       null x => outputForm(1$Integer)$EX
       reduce(_* ,[l :: EX for l in x])$List(EX)

     retractable? x == 
       null x => false
       null x.rest

     retract x == 
        not one?(#x) => error "cannot convert to Lyndon word"
        x.first

     retractIfCan x ==
        retractable? x => x.first
        "failed"
      
     length x ==
        n: Integer := +/[ length l for l in x]
        n::NNI

     recursif(x, y) ==
       null y => false
       null x => true
       x.first = y.first => recursif(rest(x), rest(y))
       lexico(x.first, y.first)

     x < y == 
       lx: NNI := length x; ly: NNI := length y 
       lx = ly => recursif(x,y)
       lx < ly

import OrderedSet
import CommutativeRing
)abbrev domain XPBWPOLY XPBWPolynomial
++ Author: Michel Petitot (petitot@lifl.fr).
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This domain constructor implements polynomials in non-commutative
++ variables written in the Poincare-Birkhoff-Witt basis from the
++ Lyndon basis.
++ These polynomials can be used to compute Baker-Campbell-Hausdorff
++ relations. \newline Author: Michel Petitot (petitot@lifl.fr).

XPBWPolynomial(VarSet:OrderedSet,R:CommutativeRing): XDPcat == XDPdef where

  WORD   ==> OrderedFreeMonoid(VarSet)
  LWORD  ==> LyndonWord(VarSet)
  LWORDS ==> List LWORD
  BASIS  ==> PoincareBirkhoffWittLyndonBasis(VarSet)
  TERM   ==> Record(k:BASIS, c:R)
  LTERMS ==> List(TERM)
  LPOLY  ==> LiePolynomial(VarSet,R)  
  EX     ==> OutputForm
  XDPOLY ==> XDistributedPolynomial(VarSet,R)
  XRPOLY ==> XRecursivePolynomial(VarSet,R)
  TERM1  ==> Record(k:LWORD, c:R)
  NNI    ==> NonNegativeInteger
  I      ==> Integer
  RN     ==> Fraction(Integer)

  XDPcat == Join(XPolynomialsCat(VarSet,R), FreeModuleCat(R, BASIS)) with
    coerce      : LPOLY -> $
      ++ \axiom{coerce(p)} returns \axiom{p}. 
    coerce      : $ -> XDPOLY
      ++ \axiom{coerce(p)} returns \axiom{p} as a distributed polynomial. 
    coerce      : $ -> XRPOLY
      ++ \axiom{coerce(p)} returns \axiom{p} as a recursive polynomial.
    LiePolyIfCan: $ -> Union(LPOLY,"failed")
      ++ \axiom{LiePolyIfCan(p)} return  \axiom{p} if \axiom{p} is a Lie polynomial.
    product     : ($,$,NNI) -> $           -- produit tronque a l'ordre n
      ++ \axiom{product(a,b,n)} returns \axiom{a*b} (truncated up to order \axiom{n}).

    if R has Module(RN) then
       exp      : ($,NNI) -> $
          ++ \axiom{exp(p,n)} returns the exponential of \axiom{p} 
          ++ (truncated up to order \axiom{n}).
       log      : ($,NNI) -> $
          ++ \axiom{log(p,n)} returns the logarithm of \axiom{p}
          ++ (truncated up to order \axiom{n}).

  XDPdef == FreeModule1(R,BASIS) add
       import(TERM)

    -- Representation
       Rep:= LTERMS 

    -- local functions
       prod1: (BASIS, $) -> $
       prod2: ($, BASIS) -> $
       prod : (BASIS, BASIS) -> $

       prod11: (BASIS, $, NNI) -> $
       prod22: ($, BASIS, NNI) -> $

       outForm : TERM -> EX
       Dexpand : BASIS -> XDPOLY
       Rexpand : BASIS -> XRPOLY
       process : (List LWORD, LWORD, List LWORD) -> $
       mirror1 : BASIS -> $

    -- functions locales
       outForm t ==
           t.c =$R 1 => t.k :: EX
           t.k =$BASIS 1 => t.c :: EX
           t.c::EX * t.k ::EX

       prod1(b:BASIS, p:$):$ ==
         +/ [t.c * prod(b, t.k) for t in p]

       prod2(p:$, b:BASIS):$ ==
         +/ [t.c * prod(t.k, b) for t in p]
 
       prod11(b,p,n) ==
           limit: I := n -$I length b
           +/ [t.c * prod(b, t.k) for t in p| length(t.k) :: I <= limit]

       prod22(p,b,n) ==
           limit: I := n -$I length b
           +/ [t.c * prod(t.k, b) for t in p| length(t.k) :: I <= limit]

       prod(g,d) ==
         d = 1 => monom(g,1)
         g = 1 => monom(d,1)
         process(reverse ListOfTerms g, first d, rest ListOfTerms d)

       Dexpand b == 
         b = 1 => 1$XDPOLY
         */ [LiePoly(l)$LPOLY :: XDPOLY for l in ListOfTerms b]

       Rexpand b ==
         b = 1 => 1$XRPOLY
         */ [LiePoly(l)$LPOLY :: XRPOLY for l in ListOfTerms b]

       mirror1(b:BASIS):$ ==
         b = 1 => 1
         lp: LPOLY := LiePoly first b
         lp := mirror lp
         mirror1(rest b) * lp :: $

       process(gauche, x, droite) ==    -- algo du "collect process"
         null gauche => monom( cons(x, droite) pretend BASIS, 1$R)
         r1, r2 : $
         not lexico(first gauche, x) =>     -- cas facile !!!
           monom(append(reverse gauche, cons(x, droite)) pretend BASIS , 1$R)

         p: LPOLY := [first gauche , x]      -- on crochete !!!
         null droite =>
           r1 :=  +/ [t.c * process(rest gauche, t.k, droite) for t in _
                      ListOfTerms p]
           r2 :=  process( rest gauche, x, list first gauche)
           r1 + r2 
         rd: List LWORD := rest droite; fd: LWORD := first droite
         r1 := +/ [t.c * process(list t.k, fd, rd) for t in  ListOfTerms p] 
         r1 := +/ [t.c * process(rest gauche, first t.k, rest ListOfTerms(t.k))_
                  for t in  r1] 
         r2 := process([first gauche, x], fd, rd)
         r2 := +/ [t.c * process(rest gauche, first t.k, rest ListOfTerms(t.k))_
                  for t in  r2]
         r1 + r2

    -- definitions
       1 == monom(1$BASIS, 1$R)

       coerce(r:R):$ == [[1$BASIS , r]$TERM ]

       coerce(p:$):EX ==
         null p => (0$R) :: EX
         le : List EX := nil
         for rec in p repeat le := cons(outForm rec, le)
         reduce(_+, le)$List(EX)

       coerce(v: VarSet):$ == monom(v::BASIS , 1$R)
       coerce(p: LPOLY):$ ==
          [[t.k :: BASIS , t.c ]$TERM for t in ListOfTerms p]

       coerce(p:$):XDPOLY ==
         +/ [t.c * Dexpand t.k for t in p]

       coerce(p:$):XRPOLY ==
         p = 0 => 0$XRPOLY
         +/ [t.c * Rexpand t.k for t in p]

       constant? p == (null p) or (leadingMonomial(p) =$BASIS 1)
       constant p == 
         null p => 0$R
         p.last.k = 1$BASIS => p.last.c
         0$R

       quasiRegular? p == (p=0) or (p.last.k ~= 1$BASIS)
       quasiRegular p == 
         p = 0 => p
         p.last.k = 1$BASIS => delete(p, maxIndex p)
         p
    
       x:$ * y:$ ==
         y = 0$$ => 0
         +/ [t.c * prod1(t.k, y) for t in x]

--       ListOfTerms p == p pretend LTERMS

       varList p == 
          lv: List VarSet := "setUnion"/ [varList(b.k)$BASIS for b in p]
          sort(lv)

       degree(p) ==
          p=0 => error "null polynomial"
          length(leadingMonomial p)

       trunc(p, n) ==
         p = 0 => p
         degree(p) > n => trunc( reductum p , n)
         p

       product(x,y,n) ==
         x = 0 => 0
         y = 0 => 0
         +/ [t.c * prod11(t.k, y, n) for t in x]

       if R has Module(RN) then
         exp (p,n) ==
             p = 0 => 1
             not quasiRegular? p => 
               error "a proper polynomial is required"
             s : $ := 1 ; r: $ := 1                  -- resultat
             for i in 1..n repeat
                k1 :RN := 1/i
                k2 : R := k1 * 1$R
                s := k2 * product(p, s, n)
                r := r + s
             r
  
         log (p,n) ==
             p = 1 => 0
             p1: $ := 1 - p
             not quasiRegular? p1 => 
               error "constant term <> 1, impossible log "
             s : $ := - 1 ; r: $ := 0                 -- resultat
             for i in 1..n repeat
               k1 :RN := 1/i
               k2 : R := k1 * 1$R
               s := product(p1, s, n)
               r := k2 * s + r
             r
 
       LiePolyIfCan p ==
         p = 0 => 0$LPOLY
         "and"/ [retractable?(t.k)$BASIS for t in p] =>
            lt : List TERM1 := _
                 [[retract(t.k)$BASIS, t.c]$TERM1 for t in p]
            lt pretend LPOLY
         "failed"

       mirror p ==
         +/ [t.c * mirror1(t.k) for t in p]

import OrderedSet
import CommutativeRing
import Module
)abbrev domain LEXP LieExponentials
++ Author: Michel Petitot (petitot@lifl.fr).
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ Management of the Lie Group associated with a
++ free nilpotent Lie algebra. Every Lie bracket with 
++ length greater than \axiom{Order} are
++ assumed to be null.
++ The implementation inherits from the \spadtype{XPBWPolynomial}
++ domain constructor: Lyndon
++ coordinates are exponential coordinates 
++ of the second kind. \newline Author: Michel Petitot (petitot@lifl.fr).

LieExponentials(VarSet, R, Order): XDPcat == XDPdef where

  EX     ==> OutputForm
  PI     ==> PositiveInteger
  NNI    ==> NonNegativeInteger
  I      ==> Integer
  RN     ==> Fraction(I)
  R      : Join(CommutativeRing, Module RN)
  Order  : PI 
  VarSet : OrderedSet
  LWORD  ==> LyndonWord(VarSet)
  LWORDS ==> List LWORD
  BASIS  ==> PoincareBirkhoffWittLyndonBasis(VarSet)
  TERM   ==> Record(k:BASIS, c:R)
  LTERMS ==> List(TERM)
  LPOLY  ==> LiePolynomial(VarSet,R)  
  XDPOLY ==> XDistributedPolynomial(VarSet,R)
  PBWPOLY==> XPBWPolynomial(VarSet, R)
  TERM1  ==> Record(k:LWORD, c:R)
  EQ     ==> Equation(R)

  XDPcat == Group with
    exp         : LPOLY -> $
      ++ \axiom{exp(p)} returns the exponential of \axiom{p}.
    log         : $ -> LPOLY
      ++ \axiom{log(p)} returns the logarithm of \axiom{p}.
    ListOfTerms : $ -> LTERMS
      ++ \axiom{ListOfTerms(p)} returns the internal representation of \axiom{p}.
    coerce      : $ -> XDPOLY
      ++ \axiom{coerce(g)} returns the internal representation of \axiom{g}.
    coerce      : $ -> PBWPOLY
      ++ \axiom{coerce(g)} returns the internal representation of \axiom{g}.
    mirror      : $ -> $
      ++ \axiom{mirror(g)} is the mirror of the internal representation of \axiom{g}.
    varList     : $ -> List VarSet
      ++ \axiom{varList(g)} returns the list of variables of \axiom{g}. 
    LyndonBasis : List VarSet -> List LPOLY
      ++ \axiom{LyndonBasis(lv)} returns the Lyndon basis of the nilpotent free
      ++ Lie algebra.
    LyndonCoordinates: $ -> List TERM1
      ++ \axiom{LyndonCoordinates(g)} returns the exponential coordinates of \axiom{g}.
    identification: ($,$) -> List EQ
      ++ \axiom{identification(g,h)} returns the list of equations \axiom{g_i = h_i},
      ++ where \axiom{g_i} (resp. \axiom{h_i}) are exponential coordinates 
      ++ of \axiom{g} (resp. \axiom{h}). 

  XDPdef == PBWPOLY add

    -- Representation
       Rep := PBWPOLY 

    -- local functions
       compareTerm1s: (TERM1, TERM1) -> Boolean
       out: TERM1 -> EX
       ident: (List TERM1, List TERM1) -> List EQ

    -- functions locales
       ident(l1, l2) ==
         import(TERM1)
         null l1 => [equation(0$R,t.c)$EQ for t in l2]
         null l2 => [equation(t.c, 0$R)$EQ for t in l1]        
         u1 : LWORD := l1.first.k; c1 :R := l1.first.c
         u2 : LWORD := l2.first.k; c2 :R := l2.first.c
         u1 = u2 =>
            r: R := c1 - c2
            r = 0 => ident(rest l1, rest l2) 
            cons(equation(c1,c2)$EQ , ident(rest l1, rest l2))
         lexico(u1, u2)$LWORD =>
            cons(equation(0$R,c2)$EQ , ident(l1, rest l2))
         cons(equation(c1,0$R)$EQ , ident(rest l1, l2))

       -- ordre lexico decroissant
       compareTerm1s(u:TERM1, v:TERM1):Boolean == lexico(v.k, u.k)$LWORD

       out(t:TERM1):EX ==
         t.c =$R 1 => char("e")$Character :: EX ** t.k ::EX
         char("e")$Character :: EX ** (t.c::EX * t.k::EX)
 
    -- definitions
       identification(x,y) ==
          l1: List TERM1 := LyndonCoordinates x
          l2: List TERM1 := LyndonCoordinates y
          ident(l1, l2)
 
       LyndonCoordinates x ==
         lt: List TERM1 := [[l::LWORD, t.c]$TERM1 for t in ListOfTerms x | _
                             (l := retractIfCan(t.k)$BASIS) case LWORD ] 
         lt := sort(compareTerm1s,lt)

       x:$ * y:$ == product(x::Rep, y::Rep, Order::I::NNI)$Rep

       exp p == exp(p::Rep , Order::I::NNI)$Rep

       log p == LiePolyIfCan(log(p,Order::I::NNI))$Rep :: LPOLY

       coerce(p:$):EX ==
          p = 1$$ => 1$R :: EX
          lt : List TERM1 := LyndonCoordinates p 
          reduce(_*, [out t for t in lt])$List(EX)


       LyndonBasis(lv) == 
         [LiePoly(l)$LPOLY for l in LyndonWordsList(lv,Order)$LWORD]

       coerce(p:$):PBWPOLY == p::Rep

       inv x ==
         x = 1 => 1
         lt:LTERMS := ListOfTerms mirror x
         lt:= [[t.k, (odd? length(t.k)$BASIS => - t.c; t.c)]$TERM for t in lt ]
         lt pretend $

