File r37/packages/ncpoly/ncdip.red artifact f5699af377 part of check-in 3c4d7b69af


module ncdip; % Non-commutative distributive polynomials.

% Author: H. Melenk, ZIB Berlin, J. Apel, University of Leipzig.


fluid '(
        ncdipvars!*         % variable set

        ncdipbase!*         % vector:
                            % the i-th entry is a list (j1,j2...)
			    % where j1,j2 ... < i
			    %    and x_i * x_j neq x_j*x_i

        ncdiptable!*        % 2-dim array:
			    % then entry (i,j) keeps the powers of the
			    % commutator [x_i,x_j] where j<i

        ncdipcircular!*     % t if one variable appears in more than one
                            % commutator

        vdpsortmode!* 
        !*gsugar dipvars!*);

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;


REDUCE Historical
REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]