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;