module decompos; % Decomposition of polynomials f(x) = g(h(x)).
% Author: Herbert Melenk <melenk@sc.zib-berlin.de>.
% Algorithms: 1. univariate case:
% V.S. Alagar, M.Tanh: Fast Polynomial Decomposition
% Algorithms, EUROCAL 1985, pp 150-153 (Springer).
%
% 2. multivariate lifting:
% J. von zur Gathen: Functional Decomposition of Polynomials:
% the Tame Case, J. Symbolic Computation (1990) 9, 281-299.
% Copyright (c) 1990 ZIB.
%
% 1-July-93 Replaced gensym calls by local name generator.
% Otherwise decompose may produce different results
% for identical input.
% 29-Apr.-93: completed normalization of multivariate results:
% shifting sign and content (field: leading coefficient)
% and absolute term to the 1st form.
global '(decomposegensym!*);
put('decompose,'psopfn,'decomposesf);
symbolic procedure decomposesf f;
'list . reverse decomposef2(simp reval car f,t)
where !*factor=nil,!*exp=t;
symbolic procedure decomposef1(f,msg);
decomposef2(f ./ 1 ,msg);
symbolic procedure decomposef2(f,msg);
begin scalar hvars,r,rr,x,y,u,vars,newvars,d;
decomposegensym!*:=1000;
vars := decomposesfvars(numr f,nil);
newvars := for each x in vars collect decomposegensym();
d := denr f;
if not domainp d
then rerror(poly,18,typerr(prepsq f,"polynomial"));
f := numr subf(numr f,pair(vars,newvars));
if length vars = 1 then r := decomposesfuni0 f
else r := decomposesfmulti(f,newvars);
hvars := '(u v w a b c d e);
for each x in vars do hvars := delete (x,hvars);
while r do
<<if cdr r
then <<y := x; x := nil;
while null x do
if hvars then <<x := car hvars; hvars := cdr hvars;
if not(x=reval x) then x := nil>>
else x:=decomposegensym();
u := prepsq subsq(car r,list(mvar numr car r . x));
if d neq 1 then<<u:=list('QUOTIENT,u,prepf d);d:=1>>;
rr := (if y then list('EQUAL,y,u) else u) . rr>>
else <<u := prepsq car r;
y := x;
rr := (if y then list('EQUAL,y,u) else u) . rr>>;
r := cdr r>>;
rr := subla(pair(newvars,vars),car rr) . cdr rr;
return rr
end;
symbolic procedure decomposesfvars(f,v);
% Select the kernels from a standard form.
if domainp f then v else
decomposesfvars(red f,
decomposesfvars(lc f,
if not member(mvar f,v)
then append(v,list mvar f) else v));
symbolic procedure decomposesfuni0 f;
for each p in decomposesfuni f collect (p ./ 1);
symbolic procedure decomposesfuni f;
% Univariate variant.
begin scalar x,y,res,ddfl,h,testf;
integer n;
n := ldeg f;
if primep n then return list f;
x := mvar f; y := decomposegensym();
ddfl := decomposefctrf decomposedf(f,x);
if length ddfl > 1 then
for each d in ddfl do
if null res and 0=remainder(n , (ldeg d + 1)) then
<<h := numr decomposeint(d,x);
if null testf then
testf := addf(f,negf numr subf(f,list(x . y)));
if quotf (testf,
addf(h,negf numr subf(h,list(x . y)))) then
res := list(decomposebacksubstuni(f,h,x),h);
if res and ldeg car res<2 then res:=nil;
>>;
if null res then return list f else
return for each u in res join decomposesfuni u
end;
symbolic procedure decomposefctrf f;
% Generate all factors of f by combining the prime factors.
begin scalar u,w,q;
q := fctrf f; u:= cdr q;
if length u = 1 and cdar u=1 then return list f;
% eliminate the two trivial factors.
w := delete(quotf(f,car q),decomposefctrf1 u);
w := delete(1,w);
return w;
end;
symbolic procedure decomposefctrf1 v;
% Collect all possible crossproducts from v.
if null v then '(1) else
begin scalar r,c,q;
c:=car v;
r:=decomposefctrf1 cdr v;
q:=for i:=1:cdr c collect exptf(car c,i);
return
append(r,
for each u in q join
for each p in r collect
multf(u,p) );
end;
symbolic procedure decomposebacksubstuni(f,h,x);
begin scalar c,g,n,p,pars,ansatz,eqs;
p := 1; n := ldeg f/ldeg h;
for i:=0:n do
<<c := mkid('coeff,i);
pars := c . pars;
ansatz := addf(multf(numr simp c,p) , ansatz);
p := multf(p,h);
>>;
pars := reverse pars;
ansatz := addf(f , negf ansatz);
eqs := decomposecoeff(ansatz,list x);
eqs := solveeval list('list . for each u in eqs collect prepf u,
'list . pars);
eqs := cdr cadr eqs; % select the only solution.
for i:= 0:n do
g := addf(g,numr simp list('times,list('expt,x,i),
caddr nth(eqs,i+1)));
return g
end;
symbolic procedure decomposedf(f,x);
% Differentiate a polynomial wrt top-level variable x.
% Returns a standard form.
if domainp f or not(mvar f = x) then nil else
if ldeg f = 1 then lc f else
mvar f .** (ldeg f - 1) .* multf(lc f,ldeg f)
.+ decomposedf(red f,x);
symbolic procedure decomposeint(f,x);
% Integrate a polynomial (standard form) wrt the (main-)variable x.
% Returns a standard quotient.
if null f then nil ./ 1 else
if domainp f then (x .** 1 .* f .+ nil) ./ 1 else
addsq(multsq((x .** (ldeg f + 1) .* 1 .+ nil)./ 1 ,
multsq(lc f./1,1 ./ldeg f+1))
, decomposeint(red f,x));
symbolic procedure decomposecoeff(f,vars);
% Select the coefficients of f wrt vars.
begin scalar o;
o := setkorder vars;
f := reorder f;
setkorder o;
return decomposecoeff1(f,vars)
end;
symbolic procedure decomposecoeff1(f,vars);
if domainp f then nil else
if not member(mvar f,vars) then list f else
nconc(decomposecoeff1(lc f,vars),decomposecoeff1(red f,vars));
symbolic procedure decomposetdg f;
% calculate total degree
if domainp f then 0 else
max(ldeg f + decomposetdg lc f, decomposetdg red f);
symbolic procedure decomposedegr(f,vl);
if domainp f then vl else
<<if ldeg f > cdr v then cdr v := ldeg f;
decomposedegr(lc f,vl);
decomposedegr(red f,vl);
vl>> where v = assoc(mvar f,vl);
symbolic procedure compose (u,v);
% Calculate f(x)=u(v(x)) for standard forms u,v.
if domainp u then u else
numr subf(u,list(mvar u . prepf v));
% Multivariate polynomial decomposition.
%
% Technique:
% select a field as domain (rational),
% map f to a strongly monic polynomial by variable transform,
% map f to a univariate image,
% decompose the univariate polynomial,
% lift decomposition to multivariate,
% convert back to original variables,
% transform back to original domain (if possible).
symbolic procedure decomposesfmulti(f,vars);
% Multivariant case: map to field (rationals).
begin scalar dm,ft,r,rr,a,q,c,p1,p2;
if null dmode!* or not flagp(dmode!*,'field) then
<<setdmode('rational,t) where !*msg=nil; dm := t;
ft := !*q2f resimp !*f2q f>> else ft := f;
r := decomposesfmulti1(ft,vars);
if dm then setdmode('rational,nil) where !*msg=nil;
if null cdr r then return list(f./1);
% if null dm then return
% for each p in r collect (p ./ 1);
% Convert back to integer polynomials.
rr := for each p in reverse r collect simp prepf p;
r := nil;
while rr and cdr rr do
<<p1 := car rr; p2 := cadr rr;
% Propagate absolute term and content from p1 to p2.
q := denr p1; a := numr p1;
while not domainp a do a := red a;
p1 := addf(numr p1,negf a);
c := decomposenormfac p1;
p1 := multsq(p1 ./ 1, 1 ./ c);
p2 := subsq(p2,list(mvar numr p2 .
list('quotient,
list('plus,list('times,decomposegensym(),prepf c),
prepf a),
prepf q)));
r := p1 . r; rr := p2 . cddr rr>>;
return car rr . r;
end;
symbolic procedure decomposesfmulti1(f,vars);
% Multivariate case: map to strongly monic polynomial.
begin scalar lvars,ft,rt,x1,a0,kord,u,sigma;
integer n,m;
% get the variable with highest degree as main variable.
u := decomposedegr(f,for each x in vars collect (x. 0));
n := -1;
for each x in u do
if n<cdr x then <<n:=cdr x; x1 := car x>>;
if n<2 then return list f;
vars := x1 . delete(x1,vars);
kord := setkorder vars;
f := reorder f;
% Convert f to a strongly monic polynomial.
n := decomposetdg f;
x1 := car vars;
lvars := for each x in cdr vars collect (x . decomposegensym());
again:
if m>10 then << rt := list f; goto ret>>;
% construct transformation sigma
sigma := for each x in lvars collect x . random 1000;
ft := numr subf(f,for each x in sigma collect
(caar x . list('plus,cdar x,list('times,x1,cdr x))));
if not domainp lc ft then <<m:=m+1; goto again>>;
a0 := lc ft; ft := quotf(ft,a0);
rt := decomposesfmnorm(ft,n,sublis(lvars,vars));
if cdr rt then
% Transform result back.
<<rt := reverse rt;
rt := numr subf(car rt,for each x in sigma collect
(cdar x . list('difference,caar x,list('times,cdr x,x1))))
. multf(a0,cadr rt) . cddr rt;
>> else rt := list f;
ret:
setkorder kord;
rt := for each p in rt collect reorder p;
% try further decomposition of central polynomial.
return if cdr rt and decomposetdg car rt>1 then
append(reverse cdr rt,decomposesfmulti1(car rt,vars))
else reverse rt;
end;
symbolic procedure decomposelmon f;
% Extract the variables of the leading monomial.
if domainp f then nil else
mvar f . decomposelmon lc f;
symbolic procedure decomposenormfac p1;
if null dmode!* or not flagp(dmode!*,'field) then
multf(numr mkabsfd decomposecont p1,decomposesign p1)
else <<while not domainp p1 do p1:=lc p1; p1>>;
symbolic procedure decomposecont f;
% Calculate the content of f if the domain is a ring.
if domainp f then f else
gcdf(decomposecont lc f, decomposecont red f);
symbolic procedure decomposesign f;
% Compute a unit factor c such that the leading coefficient of
% f/c is a positive integer.
if domainp f then numr quotsq(f ./ 1,mkabsfd f)
else decomposesign lc f;
symbolic procedure decomposesfmnorm(f,n,vars);
% Multivariate case: map strongly monic polynomial to univariate
% and lift result.
begin scalar x,x1,f0,g,u,abort,h,k,tt,q,v;
integer r,s;
x1 := car vars;
% Step 1.
f0 := numr subf(f,for each y in cdr vars collect (y . 0));
u := decomposesfuni f0;
% For multivariate we accept degree=1 polynomials as nontrivial
% but inhibit recursion.
if null cdr u then <<u:=append(u,list !*k2f x1)>>;
x := decomposegensym();
g := numr subf(car u,list (x1 . x));
r := ldeg g;
h := cadr u; u := cddr u;
while u do
<<v := car u; u:= cdr u; h := numr subf(h,list(x1 . x));
h := compose(h,v); >>;
% Step 2.
s := divide(n,r);
if not(cdr s=0) then goto fail else s := car s;
k := h;
tt := compose(decomposedf(g,x),h);
% Step 3: Hensel lifting in degree steps.
for i:=1:s do
if not abort then
% Step 4: loop step.
<<u := decomposehomog(addf(f,negf compose(g,k)),x1,i);
q := quotf(u,tt);
if u and null q then abort:=t else<<h:=q; k:=addf(k,h)>>
>>;
if abort then goto fail;
% Step 5: test result and loop for lower part.
h := k;
if f = compose(g,h) then return list(g,h);
fail: % Exit: no decomposition found.
return list f;
end;
symbolic procedure decomposehomog(f,x,d);
% F is a polynomial (standard form) in x and some other
% variables. Select that part of f, where the coefficients
% of x are monomials in total degree d.
% Result is the sum (standard form) of these monomials.
begin scalar u,v;
u := decomposehomog1(f,x,d);
for each m in u do v := addf(v,m);
return v;
end;
symbolic procedure decomposehomog1(f,x,d);
% Select the monomials.
if d<0 or null f then nil else
if domainp f then (if d=0 then list f else nil)
else begin scalar u1,u2;
u1:= decomposehomog1(lc f,x,if mvar f = x then d
else d-ldeg f);
u2:= decomposehomog1(red f,x,d);
return
nconc(
for each v in u1 collect
multf(mvar f .** ldeg f .*1 .+ nil , v),
u2);
end;
symbolic procedure decomposegensym();
compress(append('(!! !D !! !c !! !.),
explode2(decomposegensym!*:=decomposegensym!*+1)));
endmodule;
end;