Artifact 27588e8fc03a0acad32db300c4d30ac2cbc52d4c840139a0ea2e4429d0254384:
- Executable file
r38/packages/ncpoly/ncfactor.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 13626) [annotate] [blame] [check-ins using] [more...]
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;