Check-ins Using Download Hex

Artifact 6af3f8779fa13dfaac4c176b07eb77bbc6a59e893558f24fd79019c215b8ca2f:

module ncfactor;  % factorization for non-commutative polynomials.

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

fluid '(nc_factor_time nc_factor_time!* !*trnc !*ncg!-right 
        !*bcsubs2 !*gsugar ncpi!-names!* ncmp!* !*complex vdpvars!*);

% version 1.4: using the commutative factorizer as preprocessor.

switch trnc;

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);
     goto 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 goto 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 goto 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; goto 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 goto 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 goto 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:
     %###### das muss aber die Leitvariable sein!!!
   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);
   >>;

collect:
   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 goto collect else goto 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 goto collect else goto 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};
     goto 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);>>;
   w:=(cdr solveeval{'list.s,'list.vl} where dipvars!*=nil);
     % select appropiate solution.
 loop:
   nc_factorize_timecheck();
   if null w then goto 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 occuring 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 goto loop;
   q:=assoc(car vl,soa);
   if null q or cdr q=0 then goto loop;
   sol:=soa.sol;
   if all then goto 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(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-monomials 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;
REDUCE Historical
[ REDUCE Sourceforge Project | Historical SVN Repository | GitHub Mirror | SourceHut Mirror | NotABug Mirror | Chisel Mirror | Chisel RSS ]