\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra xlpoly.spad} \author{Michel Petitot} \maketitle \begin{abstract} \end{abstract} \tableofcontents \eject \section{domain FMAGMA FreeMagma} <>= import OrderedSet import RetractableTo )abbrev domain FMAGMA FreeMagma ++ 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). FreeMagma(VarSet:OrderedSet):Public == Private where WORD ==> OrderedFreeMonoid(VarSet) EX ==> OutputForm Public == Join(OrderedSet,RetractableTo VarSet,CoercibleTo WORD) with * : ($,$) -> $ ++ \axiom{x*y} returns the tree \axiom{[x,y]}. 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?}{FreeMagma}(\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?}{FreeMagma}(\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?}{FreeMagma}(\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?}{FreeMagma}(\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$FreeMagma: 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 @ \section{domain LWORD LyndonWord} A function $f \epsilon \lbrace 0,1 \rbrace$ is called acyclic if $C(F)$ consists of $n$ different objects. The canonical representative of the orbit of an acyclic function is usually called a Lyndon Word \cite{1}. If $f$ is acyclic, then all elements in the orbit $C(f)$ are acyclic as well, and we call $C(f)$ an acyclic orbit. <>= import OrderedSet import RetractableTo import Boolean import FreeMagma )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{FreeMagma} 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,CoercibleTo OFMON,CoercibleTo FreeMagma 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}. 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 == FreeMagma(VarSet) add -- Representation Rep:= FreeMagma(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:$):FreeMagma 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] @ \section{category LIECAT LieAlgebra} <>= 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 --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 @ \section{category FLALG FreeLieAlgebra} <>= 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),CoercibleTo XDPOLY,CoercibleTo XRPOLY) 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. 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}. @ \section{package XEXPPKG XExponentialPackage} <>= 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) @ \section{domain LPOLY LiePolynomial} <>= 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 FMAGMA ==> FreeMagma(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::FMAGMA * v::FMAGMA) 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) @ \section{domain PBWLB PoincareBirkhoffWittLyndonBasis} <>= 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,CoercibleTo WORD) with 1: constant -> % ++ \spad{1} returns the empty list. 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 @ \section{domain XPBWPOLY XPBWPolynomial} <>= 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),CoercibleTo XDPOLY,CoercibleTo XRPOLY) with coerce : LPOLY -> $ ++ \axiom{coerce(p)} returns \axiom{p}. 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] @ \section{domain LEXP LieExponentials} <>= 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 == Join(Group,CoercibleTo XDPOLY,CoercibleTo PBWPOLY) 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}. 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 $ @ \section{License} <>= --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. @ <<*>>= <> <> <> <> <> <> <> <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} {\bf http://www.mathe2.uni-bayreuth.de/frib/html/canonsgif/canons.html} \end{thebibliography} \end{document}