module ncdip; % Non-commutative distributive polynomials.
% Author: H. Melenk, ZIB Berlin, J. Apel, University of Leipzig.
symbolic procedure ncdsetup!* u;
% U is a list of algebraic arguments:
% 1. list of variables,
% 2. list of commutator relations in explicit form x*y=y*x + r
% where ord(r) < ord(x*y) .
% All variable pairs whitch do not occur here are considered
% communtative.
begin scalar x,y,w,vars,lh,z,lv,r,q;
vars := for each x in cdr listeval(car u,nil)
collect reval x;
ncdipcircular!*:=nil;
if null vdpsortmode!* then vdpsortmode!*:= 'lex;
vdpinit2 (ncdipvars!*:=vars);
lv:=length vars;
ncdipbase!*:=mkvect lv;
ncdiptable!*:=mkvect lv;
for i:=1:lv do putv(ncdiptable!*,i,mkvect lv);
q:=cdr listeval(cadr u,nil);
while q do
<<r:=car q;q:=cdr q;
if (eqcar(r,'equal) or eqcar(r,'replaceby)) and
eqcar(lh:=cadr r,'times) and
null cdddr r and
(w:=member(y:=caddr lh,vars)) and
member(x:=cadr lh,w)
then
<< % does any variable other than x or y appear in the rhs?
if smember(x,q) or smember(y,q) then ncdipcircular!*:=t;
for each v in vars do
if v neq x and v neq y and smember(v,caddr r) then
ncdipcircular!*:=t;
% establish the commutator in DIP form.
w:=ncdipndx(x,vars,1).ncdipndx(y,vars,1);
r:=a2dip(z:=reval caddr r);
if evcomp(dipevlmon r,dipevlmon a2dip {'times,x,y})<0 then
typerr({'equal,{'times,x,y},z},
"commutator under current term order");
getv(ncdipbase!*,car w):= sort(cdr w.getv(ncdipbase!*,car w),'lessp);
getv(getv(ncdiptable!*,car w),cdr w):= {'(1 . 1).r}>>
else typerr(r,"commutator ")>>end;
symbolic procedure ncdipndx(x,vars,n);
if null vars then 0 else if x=car vars then n else ncdipndx(x,cdr vars,n #+ 1);
%------------ noncom multiply ----------------------------
symbolic procedure vdp!-nc!-m!*p(bc,ev,p);
% multiply polynomial p left by monomial (bc,ev).
begin scalar r,s;
r:=dip2vdp dip!-nc!-m!*p(bc,ev,vdppoly p);
if !*gsugar then
<<s:=gsugar p; for each e in ev do s:=s+e;r:=gsetsugar(r,s)>>;
return r end;
symbolic procedure vdp!-nc!-prod(u,v);
% non-commutative product of two distributive polynomials.
begin scalar r;
r:=dip2vdp dip!-nc!-prod(vdppoly u,vdppoly v);
if !*gsugar then r:=gsetsugar(r,gsugar u + gsugar v);
return r end;
symbolic procedure dip!-nc!-prod(u,v);
% We distribute first over the shorter of the two factors.
if length u < length v then dip!-nc!-prod!-distleft(u,v)
else dip!-nc!-prod!-distright(u,v);
symbolic procedure dip!-nc!-prod!-distleft(u,v);
if dipzero!? u then u else
dipsum(dip!-nc!-m!*p!-distleft(diplbc u,dipevlmon u,v),
dip!-nc!-prod!-distleft(dipmred u,v));
symbolic procedure dip!-nc!-m!*p!-distleft(bc,ev,p);
if dipzero!? p then nil else
begin scalar lev,lbc,q;
lev:=dipevlmon p;lbc:=diplbc p;
p:=dip!-nc!-m!*p!-distleft(bc,ev,dipmred p);
q:=dip!-nc!-ev!-prod(bc,ev,lbc,lev);
return dipsum(p,q)end;
symbolic procedure dip!-nc!-prod!-distright(u,v);
if dipzero!? v then v else
dipsum(dip!-nc!-m!*p!-distright(u,diplbc v,dipevlmon v),
dip!-nc!-prod!-distright(u,dipmred v));
symbolic procedure dip!-nc!-m!*p!-distright(p,bc,ev);
if dipzero!? p then nil else
begin scalar lev,lbc,q;
lev:=dipevlmon p;lbc:=diplbc p;
p:=dip!-nc!-m!*p!-distright(dipmred p,bc,ev);
q:=dip!-nc!-ev!-prod(lbc,lev,bc,ev);
return dipsum(p,q)end;
symbolic procedure dip!-nc!-ev!-prod(bc1,ev1,bc2,ev2);
% compute (bc1*ev1) * (bc2*ev2). Result is a dip.
dip!-nc!-ev!-prod1(ev1,1,dipfmon(bcprod(bc1,bc2),ev2));
symbolic procedure dip!-nc!-ev!-prod1(ev,n,r);
% loop over ev and n (counter). NOTE: ev must be processed from right to left!
if null ev then r else
dip!-nc!-ev!-prod2(car ev,n,dip!-nc!-ev!-prod1(cdr ev,n#+1,r));
symbolic procedure dip!-nc!-ev!-prod2(j,n,r);
% muliply x_n^j * r
if j=0 or dipzero!? r then r else
begin scalar ev,bc,r0,w,s,evl,evr;integer i;
ev:=dipevlmon r;bc:=diplbc r;
r:=dip!-nc!-ev!-prod2(j,n,dipmred r);
% collect the variables in ev which do not commute with x_n;
w:=getv(ncdipbase!*,n);
while w and nth(ev,car w)=0 do w:=cdr w;
% no commutator?
if null w then
<<r0:=dipfmon(bc,dipevadd1var(j,n,ev));return dipsum(r0,r)>>;
% We handle now the leftmost commutator and
% push the rest of the problem down to the recursion:
% split the monmial into parts left and
% right of the noncom variable and multiply these.
w:=car w;s:=nth(ev,w);
% Split the ev into left and right part.
i:=0;for each e in ev do <<i:=i#+1;if i<w then<<evl:=e.evl;evr:=0 .evr>>
else if i=w then <<evl:=0 .evl;evr:=0 .evr>>
else <<evl:=0 .evl;evr:=e .evr>> >>;
evl:=reversip evl;evr:=reversip evr;
r0:=dip!-nc!-get!-commutator(n,j,w,s);
% multiply by left exponent
r0:=if ncdipcircular!* then
<<r0:=dip!-nc!-prod(dipfmon(a2bc 1,evl),r0);
r0:=dip!-nc!-prod(r0,dipfmon(bc,evr))>> else
<<r0:=dipprod(dipfmon(a2bc 1,evl),r0);
r0:=dipprod(r0,dipfmon(bc,evr))>>;
done:return dipsum(r0,r)end;
symbolic procedure dip!-nc!-get!-commutator(n1,e1,n2,e2);
% Compute the commutator for y^e1*x^e2 where y is
% the n1-th variable and x is the n2-th variable.
% The commutators for power products are computed
% recursively when needed. They are stored in a data base.
% I assume here that the commutator for (1,1) has been
% put into the data base before. We update the table
% in place by nconc-ing new pairs to its end.
begin scalar w,r,p;
w:=getv(getv(ncdiptable!*,n1),n2);p:=e1.e2;
if (r:=assoc(p,w)) then return cdr r;
% compute new commutator recursively:
% first e1 downwards, then e2
r:=if e1>1 then
% compute y^e1*x^e2 as y*(y^(e1-1)*x^e2)
dip!-nc!-ev!-prod2(1,n1,dip!-nc!-get!-commutator(n1,e1#-1,n2,e2))
else
% compute y*x^e2 as (y*x^(e2-1))*x
dip!-nc!-prod(dip!-nc!-get!-commutator(n1,1,n2,e2#-1),dipfvarindex n2);
nconc(w,{(p.r)});
return r end;
symbolic procedure dipfvarindex n;
% Make a dip from a single variable index.
a2dip nth(dipvars!*,n);
symbolic procedure dipevadd1var(e,n,ev);
% add e into the nth position of ev.
if null ev or n<1 then ev else
if n=1 then (car ev #+ e).cdr ev else car ev.dipevadd1var(e,n#-1,cdr ev);
% ------------ conversion algebraic => nc dip --------------
symbolic procedure a2ncdip a;
if atom a then a2dip a else
if car a = 'difference then
a2ncdip{'plus,cadr a,{'times,-1,caddr a}} else
if car a = 'minus then
a2ncdip{'times,-1,cadr a} else
if car a='expt and fixp caddr a then
a2ncdip('times.for i:=1:caddr a collect cadr a) else
if car a='plus then
begin scalar r;
r:=a2ncdip cadr a;
for each w in cddr a do r:=dipsum(r,a2ncdip w);
return r end else
if car a='times then
begin scalar r;
r:=a2ncdip cadr a;
for each w in cddr a do r:=dip!-nc!-prod(r,a2ncdip w);
return r end else
if car a='!*sq then a2ncdip prepsq cadr a else a2dip a;
symbolic procedure a2ncvdp a;dip2vdp a2ncdip a;
endmodule;;end;