reduce-historical/r38/packages/ncpoly/ncfactor.red

module ncfactor;  % factorization for non-commutative polynomials.

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

% version 1.4: using the commutative factorizer as preprocessor.
% Oct 2001: using "sove", hoping, that the user did not switch off 'varopt'.

share nc_factor_time;    % time limit in milliseconds.

nc_factor_time:=0;

algebraic operator cc!*;

symbolic procedure nc_factorize u;
 begin scalar r,o,!*gsugar,comm,cr,cl;
   o:=apply1('torder,'(gradlex));
   nc!-gsetup();
   comm := nc_commfactors!* u;
   cl:=car comm; u:=cadr comm;cr:=caddr comm; 
   if constant_exprp u then (if u neq 1 then cl:=u.cl)
      else
     r:=for each p in nc_factorize0(a2ncvdp u,nil,nil,nil,nil,nil)
          collect num vdp2a p;
   o:=apply1('torder,{o});
   return'list.append(cl,append(r,cr))end;

symbolic operator nc_factorize;

% copyd('nc_commfactors!*,'nc_commfactors);
symbolic procedure nc_commfactors u;
 begin scalar o,!*gsugar,comm,cr,cl;
   o:=apply1('torder,'(gradlex));
   nc!-gsetup();
   comm:=nc_commfactors!* u;
   cl:=car comm;u:=cadr comm;cr:=caddr comm; 
   o:=apply1('torder,{o});
   return{'list,'list.cl,u,'list.cr}end;

symbolic operator nc_commfactors;

symbolic procedure nc_commfactors!* u;
 (begin scalar f,ff,uu,comm,l,crl,cll,!*ncg!-right,w; 
   uu:=sublis(ncpi!-names!*,numr simp u);
   comm:=(fctrf reorder uu) where ncmp!*=nil;
   if null cddr comm and cdadr comm=1 then 
   <<if !*trnc then writepri("no commutative factors found",'only);
     go to no_comm>>;
   l:=for each f in cdr comm join
     for i:=1:cdr f collect reval prepf car f;
   if !*trnc then writepri("testing commutative factors:",'only);
   uu:=a2ncvdp u;
   while l do
   <<f:=car l;l:=cdr l;
     if !*trnc then writepri(mkquote f,'first); 
     !*ncg!-right:=right;
     if vdpzero!? cdr(w:=nc!-qremf(uu,ff:=a2ncvdp f))then
     <<if !*trnc then writepri(nc_dir(),'last);cll:=append(cll,{f});uu:=car w>>
     else
     if vdpzero!? cdr<<!*ncg!-right:=not right;w:=nc!-qremf(uu,ff)>>
       then<<if !*trnc then writepri(nc_dir(),'last);crl:=f.crl;uu:=car w>>
     else if !*trnc then writepri(" -- discarded",'last)>>;  
  if null crl and null cll then go to no_comm;  
  u:=vdp2a uu; 
  if !*trnc then 
    <<writepri("remaining noncom  part:",'first);writepri(mkquote u,'last)>>;
 no_comm:return {crl,u,cll};
 end)where right=!*ncg!-right;

symbolic procedure nc_dir();if !*ncg!-right then " right" else " left";

symbolic procedure oneside!-factor(w,m,all);
 % NOTE: we must perform a factorization based on left
 % division (m='l) for obtaining a right factor.
 begin scalar u,d,r,mx,o,!*gsugar;
    % preprocessing for psopfn.
  d:=r:=0;
  u:=reval car w;
  if cdr w then<<d:=reval car(w:=cdr w);if cdr w then r:=reval cadr w>>;
    % preparing for the altorithm.
  o:=apply1('torder,'(gradlex));
  nc!-gsetup();
  if r=0 or r='(list)then r:=nil else 
  <<r:=cdr listeval(r,nil);
    r:=vdpevlmon a2vdp(if null cdr r then reval car r else
      'times.for each y in r collect reval y)>>; 
  d:=reval d; 
  if d=0 then d:=1000 else if not fixp d then<<mx:=vdpevlmon a2vdp d;d:=1000>>;
  r:=nc_factorize0(a2ncvdp u,m,d,r,mx,all);
  o:=apply1('torder,{o});
  return for each w in r collect num vdp2a w end;

put('left_factor,'psopfn,
     function (lambda(w);<<w:=oneside!-factor(w,'r,nil) or w;reval car w>>));

put('left_factors,'psopfn,
     function (lambda(w);'list. oneside!-factor(w,'r,t)));

put('right_factor,'psopfn,
     function (lambda(w);<<w:=oneside!-factor(w,'l,nil) or w;reval car w>>));

put('right_factors,'psopfn,
     function (lambda(w);'list.oneside!-factor(w,'l,t)));

algebraic procedure nc_factorize_all u;
  % Compute all possible factorizations based on successive
  % right factor extraction.
 begin scalar !*ncg!-right,d,f,w,wn,q,r,trnc,nc_factor_time!*;
   nc_factor_time!*:=lisp time();
   trnc:=lisp !*trnc;lisp(!*trnc:=nil);
   w:={{u}};r:={};lisp(!*ncg!-right:=nil);
 loop:if w={} then go to done;
   lisp(wn:='(list));
   for each c in w do
   <<lisp(q:= cadr c);
     f:=right_factors(q,{},{});
     if trnc then write "ncfctrall: Right factors of (",q,"): ",f;
     if f={} then r:=c.r;
     for each fc in f do
     <<d:=nc_divide(q,fc);
       if trnc then write "ncfctrall: Quotient (",q,") / (",fc,"): ",d;
       wn:=(first d.fc.rest c).wn>>>>;
   w:=wn; go to loop;
 done:lisp(!*trnc:=trnc);
   return r end;

symbolic procedure nc_factorize0(u,m,d,rs,mx,all);
 <<if not numberp nc_factor_time!* then nc_factor_time!*:=time();
   nc_factorize1(u,m,d,rs,mx,all)>>where nc_factor_time!*=nc_factor_time!*;
   
symbolic procedure nc_factorize1(u,m,d,rs,mx,all);
 % split all left(right) factor of u off.
 % u:  polynomial,
 % m:  mode: restriction for left or right factor:
 % d:  maximum degree restriction,
 % r:  variable set restriction (r is an exponent vector).
 % mx: maximum exponent for each variable (is an exponent vector).
 % all: true if we look for all right(left) factors.
 begin scalar ev,evl,evlx,f,ff,!*ncg!-right;
  nc_factorize_timecheck();
  mx:=if null mx then for each y in vdpvars!* collect 1000 else
    for each y in mx collect if y>0 then y else 1000;
  if !*trnc then<<prin2 "factorize ";vdpprint u>>;
  ev:=vdpevlmon u;
  if vevzero!? ev then return{u};
  d:=d or vevtdeg ev/2;
  evlx:=sort(nc_factorize1!-evl ev,function(lambda(x,y);vevcomp(x,y)<0));
  if m='r then go to r;
    % factors up to n
  evl:=evlx;
  while (null f or all) and evl and vevtdeg car evl<=d do
  <<if not vevzero!? car evl 
      and car evl neq ev
           % testing support;
      and(null rs or vevmtest!?(car evl,rs)) 
           % testing maximal degrees;
      and vevmtest!?(mx,car evl)
    then f:=append(f,nc_factorize2(u,car evl,rs,mx,all));
   evl:=cdr evl>>;
  if f or m='l then go to c;
    % right factors up to tdg-n
  d:=vevtdeg ev -d; 
r:!*ncg!-right:=t;
  evl:=evlx;
  while (null f or all)and evl and vevtdeg car evl<=d do
  <<if not vevzero!? car evl 
        and car evl neq ev
           % testing support;
        and(null rs or vevmtest!?(car evl,rs)) 
           % testing maximal degrees;
      and vevmtest!?(mx,car evl)
     then f:=append(f,nc_factorize2(u,car evl,rs,mx,all));
    evl:=cdr evl>>;
c:if null f then return if m then nil else{u};
  if all then return f;
    % only one factor wanted?
  if m then return{cdr f};
  ff:=nc_factorize1(car f,nil,nil,nil,mx,all);
  return if !*ncg!-right then append({cdr f},ff)else append(ff,{cdr f})end;
  
symbolic procedure nc_factorize1!-evl u;
  % Collect all monomials dividing u.
   if null u then'(nil) else
   (for i:=0:car u join
     for each e in w collect i.e)where w=nc_factorize1!-evl cdr u;

algebraic operator ncc!@;

symbolic procedure nc_factorize2(u,ev,rs,mx,all);
  begin scalar ar,p,q,vl,r,s,so,sol,w,y;integer n;
   scalar !*bcsubs2;
   nc_factorize_timecheck();
   p:=a2dip 0;
   if !*trnc then
   <<prin2 if !*ncg!-right then "right " else "left ";
     prin2 "Ansatz for leading term > ";
     vdpprin2 vdpfmon(a2bc 1,ev);
     prin2 " < time so far:";
     prin2 (time()-nc_factor_time!*);
     prin2t "ms">>;
     % establish formal Ansatz.
   for each e in nc_factorize2evl(ev,rs,mx) do
   <<q:={'ncc!@,n:=n+1};p:=dipsum(p,dipfmon(a2vbc q,e))>>;
   w:=p;
   while not dipzero!? w do<<vl:=bc2a diplbc w.vl;w:=dipmred w>>;
   vl:=reversip vl;
   p:=dip2vdp p;
     %  prin2 "complete Ansatz:";vdpprint p;
     % pseudo division.
   r:=nc!-normalform(u,{p},nil,nil);
   nc_factorize_timecheck();
   while not vdpzero!? r do<<s:=vbc2a vdplbc r.s;r:=vdpred r>>;
   if !*trnc then
   <<prin2t "internal equation system:";writepri(mkquote('list.s),'only)>>;
     % solve system
     % 1. look for a free variable:
     %###### but that must be the leading variable!!!
   for each v in vl do if not smember(v,s) then so:=v;
   if !*trnc and so then<<prin2"free:";prin2t so>>;
   if so then sol:={(so.1).for each v in vl collect v.0};
   if null sol or all then sol:=append(sol,nc_factsolve(s,vl,all));
   if null sol then return nil;
   if !*trnc then 
   <<prin2t "internal solutions:";
     for each s in so do
     <<for each q in s do
       <<writepri(mkquote car q,'first);
         writepri(mkquote " = ",nil);
         writepri(mkquote cdr q,'last)>>;
       prin2t "=====================================">>;
   % prin2 "check internal solution:";
   % for each e in s do writepri(mkquote aeval sublis(so,e),'only);
   >>;
coll:nc_factorize_timecheck();
   so:=car sol;sol:=cdr sol;
   y:=dip2vdp dippolish dipsubf(so,vdppoly p);
     % leading term preserved?
  % if vdpevlmon y neq vdpevlmon p then
   %  return nil;
   %  prin2 "computed factor:";vdpprint y;
   if vevzero!? vdpevlmon y then
      if not all then return nil else
      if sol then go to coll else go to done_all;
     % turn on bcsubs2 if there is an algebraic number.
   if smemq('expt,y) or smemq('sqrt,y) or smemq('root_of,y) then !*bcsubs2:=t;
   w:=nc!-qremf(u,y);
   if not vdpzero!? cdr w then
    <<prin2 "division failure";
      vdpprint u;prin2t "/";
      vdpprint y;prin2 "=> ";vdpprint car w;
      prin2 "rem: ";vdpprint cdr w;
       rederr "noncom factorize">>;
   if !*trnc then
   <<terpri();prin2 "splitting into > ";
     vdpprin2 car w;prin2t " < and";prin2 " > ";
     vdpprin2 y;prin2t " <";terpri()>>;
   ar:=y.ar;
   if all then if sol then go to coll else go to done_all;
done_one:return car w.y;
done_all:return ar end;

symbolic procedure nc_factsolve(s,vl,all);
  begin scalar v,sb,ns,so,soa,sol,nz,w,q,z,r,abort;
    % 1st phase: divide out leading term variable,
    % remove zero products, and terminate for explicitly
    % unsolvable system.
   v:=numr simp car vl;
   ns:=for each e in s collect numr simp e;
    % remove factors of leading coefficient, 
    % remove trivial parts and propagate them into system.
   r:=t;
   while r do
   <<r:=nil; s:=ns; ns:=nil;
     for each e in s do if not abort then
     <<e:=absf numr subf(e,sb);
       while(q:=quotf(e,v))do e:=q;
       if null e then nil else
       if domainp e or not(mvar e member vl)then abort:=t else
       if null red e and domainp lc e then
       <<w:=mvar e;sb:=(w.0).sb;r:=t;vl:=delete(w,vl)>>
       else if not member(e,ns)then ns:=e.ns>>>>;
   if abort or null vl then return nil;
   nc_factorize_timecheck();
     % all equations solved, free variable(s) left
   if null ns and vl then 
   <<sol:={for each x in vl collect x.1};go to done>>;
     % solve the system.
   s:=for each e in ns collect prepf e;
   if !*trnc then
    <<prin2 "solving ";
      prin2 length s;prin2 " polynomial equations for ";
      prin2 length vl;
      prin2t"variables";
      for each e in s do writepri(mkquote e,'only)>>;
        % modification HM 24.10.2001: introduction of the fluid variable
        % '*varoptt' and setting it 't' locally.
   w:=(cdr solveeval{'list.s,'list.vl} where dipvars!*=nil);
     % Select appropriate solution.
 loop:nc_factorize_timecheck();
   if null w then go to done;
   so:=cdr car w;w:=cdr w;soa:=nil;
   if smemq('i,so)and null !*complex then go to loop;
     % Insert values for non occurring variables.
   for each y in vl do if not smember(y,so)then<<soa:=(y.1).soa; nz:=t>>;
   for each y in so do
   <<z:=nc_factorize_unwrap(reval caddr y,soa); 
     nz:=nz or z neq 0;
     soa:=(cadr y.z).soa>>;
     % don't accept solution with leading term 0.
   if not nz then go to loop;
   q:=assoc(car vl,soa);
   if null q or cdr q=0 then go to loop;
   % Make sure solutions are in lowest terms.
   soa:=for each j in soa collect(car j.sublis(soa,cdr j));
   sol:=soa.sol;
   if all then go to loop;
 done:sol:=for each s in sol collect append(sb,s);
   if !*trnc then
    <<prin2t "solutions:";
      for each w in sol do
       writepri(mkquote('list.
         for each s in w collect{'equal,car s,cdr s}),'only);
      prin2t "-------------------------">>;
   return sol end;

symbolic procedure dipsubf(a,u);
  % construct polynomial u with coefficients from a.
 if dipzero!? u then nil else
  <<q:=if q then cdr q else diplbc u;
    if q neq 0 then dipmoncomp(a2bc q,dipevlmon u,r) else r>>
      where q=assoc(bc2a diplbc u,a),r=dipsubf(a,dipmred u);
        
symbolic procedure dippolish p1;diprectoint(p1,diplcm p1);
 
symbolic procedure nc_factorize_unwrap(u,s);
   if atom u then u else
   if eqcar(u,'arbcomplex)then 1 else
   (if q then cdr q else
   for each x in u collect nc_factorize_unwrap(x,s))where q=assoc(u,s);

symbolic procedure nc_factorize2evl(ev,rs,mx);
  % make list of monomials below ev in gradlex ordering,
  % but only those which occur in rs (if that is non-nil)
  % and which have the maximal degress of mx.
 for each q in nc_factorize2!-evl1(min(evtdeg mx,evtdeg ev),length ev,rs) 
  join if not vevcompless!?(ev,q) and vevmtest!?(mx,q)then{q}; 

symbolic procedure nc_factorize2!-evl1(n,m,rs);
% Collect all 'm' exponent vectors with total degree <='n'.
   if m=0 then'(nil)else
   for i:=0:(if null rs or car rs>0 then n else 0)join
     for each e in nc_factorize2!-evl1(n#-i,m#-1,if rs then cdr rs) 
       collect i.e;

symbolic procedure nc_factorize_timecheck();
   if fixp nc_factor_time and nc_factor_time>0 and
     (time() - nc_factor_time!*) > nc_factor_time
       then rederr "time overflow in noncom. factorization";

endmodule;;end;