module dpoly;
COMMENT
##################
## ##
## POLYNOMIALS ##
## ##
##################
Polynomial vectors and polynomials are handled in a unique way using
the module component of monomials to store the vector component. If
the component is 0, we have a polynomial, otherwise a vector. They
are represented in a distributive form (dpoly for short).
Informal syntax of (vector) polynomials :
<dpoly> ::= list of <term>s
<term> ::= ( <monomial> . <base coefficient> )
END COMMENT;
% ----------- constructors and selectors -------------------
symbolic procedure dp_lc p;
% Leading base coefficient of the dpoly p.
cdar p;
symbolic procedure dp_lmon p;
% Leading monomial of the dpoly p.
caar p;
symbolic procedure dp_term (a,e);
% Constitutes a term from a:base coeff. and e:monomial.
(e . a);
symbolic procedure dp_from_ei n;
% Returns e_i as dpoly.
list dp_term(bc_fi 1,mo_from_ei n);
symbolic procedure dp_fi n;
% dpoly from integer
if n=0 then nil else
list dp_term(bc_fi n,mo_zero());
symbolic procedure dp_fbc c;
% Converts the base coefficient c into a dpoly.
if bc_zero!? c then nil else
list dp_term(c,mo_zero());
% ------------ dpoly arithmetics ---------------------------
symbolic procedure dp!=comp(i,v);
if null v then nil
else if eqn(mo_comp dp_lmon v,i) then car v . dp!=comp(i,cdr v)
else dp!=comp(i,cdr v);
symbolic procedure dp_comp(i,v);
% Returns the (polynomial) component i of the vector v.
for each x in dp!=comp(i,v) collect (mo_deletecomp car x) . cdr x;
symbolic procedure dp!=mocompare (t1,t2);
% true <=> term t1 is smaller than term t2 in the current term order.
eqn(mo_compare (car t1, car t2),1);
symbolic procedure dp_neworder p;
% Returns reordered dpoly p after change of the term order.
sort(for each x in p collect (mo_neworder car x) . cdr x,
function dp!=mocompare);
symbolic procedure dp_neg p;
% Returns - p for the dpoly p.
for each x in p collect (car x . bc_neg cdr x);
symbolic procedure dp_times_mo (mo,p);
% Returns p * x^mo for the dpoly p and the monomial mo.
for each x in p collect (mo_sum(mo,car x) . cdr x);
symbolic procedure dp_times_bc (bc,p);
% Returns p * bc for the dpoly p and the base coeff. bc.
for each x in p collect (car x . bc_prod(bc,cdr x));
symbolic procedure dp_times_bcmo (bc,mo,p);
% Returns p * bc * x^mo for the dpoly p, the monomial mo and the base
% coeff. bc.
for each x in p collect (mo_sum(mo,car x) . bc_prod(bc,cdr x));
symbolic procedure dp_times_ei(i,p);
% Returns p * e_i for the dpoly p.
dp_neworder for each x in p collect (mo_times_ei(i,car x) . cdr x);
symbolic procedure dp_project(p,k);
% Delete all terms x^a*e_i with i>k.
for each x in p join if mo_comp car x <= k then {x};
symbolic procedure dp_content p;
% Returns the leading coefficient, if invertible, or the content of
% p.
if null p then bc_fi 0
else begin scalar w;
w:=dp_lc p; p:=cdr p;
while p and not bc_inv w do
<< w:=bc_gcd(w,dp_lc p); p:=cdr p >>;
return w
end;
symbolic procedure dp_mondelete(p,s);
% Returns (p.m) with common monomial factor m with support in the
% var. list s deleted.
if null p or null s then (p . mo_zero()) else
begin scalar cmf;
cmf:=dp!=cmf(p,s);
if mo_zero!? cmf then return (p . cmf)
else return
cons(for each x in p collect mo_diff(car x,cmf) . cdr x,cmf)
end;
symbolic procedure dp!=cmf(p,s);
begin scalar a;
a:=mo_seed(dp_lmon p,s); p:=cdr p;
while p and (not mo_zero!? a) do
<< a:=mo_gcd(a,mo_seed(dp_lmon p,s)); p:=cdr p >>;
return a
end;
symbolic procedure dp_unit!? p;
% Tests whether lt p of the dpoly p is a unit.
% This means : p is a unit, if the t.o. is noetherian
% or : p is a local unit, if the t.o. is a tangentcone order.
p and (mo_zero!? dp_lmon p);
symbolic procedure dp_simp pol;
% Returns (pol_new . z) with
% pol_new having leading coefficient 1 or
% dp_content pol canceled out
% and pol_old = z * dpoly_new .
if null pol then pol . bc_fi 1
else begin scalar z,z1;
if (z:=bc_inv (z1:=dp_lc pol)) then
return dp_times_bc(z,pol) . z1;
% -- now we assume that base coefficients are a gcd domain ----
z:=dp_content pol;
if bc_minus!? z1 then z:=bc_neg z;
pol:=for each x in pol collect
car x . car bc_divmod(cdr x,z);
return pol . z;
end;
symbolic procedure dp_prod(p1,p2);
% Returns p1 * p2 for the dpolys p1 and p2.
if length p1 <= length p2 then dp!=prod(p1,p2)
else dp!=prod(p2,p1);
symbolic procedure dp!=prod(p1,p2);
if null p1 or null p2 then nil
else
begin scalar v;
for each x in p1 do
v:=dp_sum( dp_times_bcmo(cdr x,car x, p2 ),v);
return v;
end;
symbolic procedure dp_sum(p1,p2);
% Returns p1 + p2 for the dpolys p1 and p2.
if null p1 then p2
else if null p2 then p1
else begin scalar sl,al;
sl := mo_compare(dp_lmon p1, dp_lmon p2);
if sl = 1 then return car p1 . dp_sum(cdr p1, p2);
if sl = -1 then return car p2 . dp_sum(p1, cdr p2);
al := bc_sum(dp_lc p1, dp_lc p2);
if bc_zero!? al then return dp_sum(cdr p1, cdr p2)
else return dp_term(al,dp_lmon p1) . dp_sum(cdr p1, cdr p2)
end;
symbolic procedure dp_diff(p1,p2);
% Returns p1 - p2 for the dpolys p1 and p2.
dp_sum(p1, dp_neg p2);
symbolic procedure dp_power(p,n);
% Returns p^n for the dpoly p.
if (not fixp n) or (n < 0) then typerr(n," exponent")
else if n=0 then dp_fi 1
else if n=1 then p
else if null cdr p then dp!=power1(p,n)
else dp!=power(p,n);
symbolic procedure dp!=power1(p,n); % For monomials.
list dp_term(bc_power(dp_lc p,n),mo_power(dp_lmon p,n));
symbolic procedure dp!=power(p,n);
if n=1 then p
else if evenp n then dp!=power(dp_prod(p,p),n/2)
else dp_prod(p,dp!=power(dp_prod(p,p),n/2));
symbolic procedure dp_tcpart p;
% Return the homogeneous degree part of p of highest degree.
if null p then nil
else begin scalar d,u; d:=car mo_deg caar p;
while p and (d=car mo_deg caar p) do
<< u:=car p . u; p:=cdr p >>;
return reversip u;
end;
symbolic procedure dp_deletecomp p;
% delete the component part from all terms.
dp_neworder for each x in p collect mo_deletecomp car x . cdr x;
symbolic procedure dp_factor p;
for each y in cdr ((fctrf numr simp dp_2a p) where !*factor=t)
collect dp_from_a prepf car y;
% ------ Converting prefix forms into dpolys ------------------
symbolic procedure dp_from_a u;
% Converts the algebraic (prefix) form u into a dpoly.
if eqcar(u,'list) or eqcar(u,'mat) then typerr(u,"dpoly")
else if atom u then dp!=a2dpatom u
else if not atom car u or not idp car u
then typerr(car u,"dpoly operator")
else (if x='dp!=fnpow then dp!=fnpow(dp_from_a cadr u,caddr u)
else if x then
apply(x,list for each y in cdr u collect dp_from_a y)
else dp!=a2dpatom u)
where x = get(car u,'dp!=fn);
symbolic procedure dp!=a2dpatom u;
% Converts the atom (or kernel) u into a dpoly.
if u=0 then nil
else if numberp u or not member(u, ring_all_names cali!=basering)
then list dp_term(bc_from_a u,mo_zero())
else list dp_term(bc_fi 1,mo_from_a u);
symbolic procedure dp!=fnsum u;
% U is a list of dpoly expressions. The result is the dpoly
% representation for the sum. Analogously for the other symbolic
% procedures below.
(<<for each y in cdr u do x := dp_sum(x,y); x>>) where x = car u;
put('plus,'dp!=fn,'dp!=fnsum);
put('plus2,'dp!=fn,'dp!=fnsum);
symbolic procedure dp!=fnprod u;
(<<for each y in cdr u do x := dp_prod(x,y); x>>) where x = car u;
put('times,'dp!=fn,'dp!=fnprod);
put('times2,'dp!=fn,'dp!=fnprod);
symbolic procedure dp!=fndif u; dp_diff(car u, cadr u);
put('difference,'dp!=fn,'dp!=fndif);
symbolic procedure dp!=fnpow(u,n); dp_power(u,n);
put('expt,'dp!=fn,'dp!=fnpow);
symbolic procedure dp!=fnneg u;
( if null v then v else dp_term(bc_neg dp_lc v,dp_lmon v) . cdr v)
where v = car u;
put('minus,'dp!=fn,'dp!=fnneg);
symbolic procedure dp!=fnquot u;
if null cadr u or not null cdadr u
or not mo_zero!? dp_lmon cadr u
then typerr(dp_2a cadr u,"distributive polynomial denominator")
else dp!=fnquot1(car u,dp_lc cadr u);
symbolic procedure dp!=fnquot1(u,v);
if null u then u
else dp_term(bc_quot(dp_lc u,v), dp_lmon u) .
dp!=fnquot1(cdr u,v);
put('quotient,'dp!=fn,'dp!=fnquot);
% -------- Converting dpolys into prefix forms -------------
% ------ Authors: R. Gebauer, A. C. Hearn, H. Kredel -------
symbolic procedure dp_2a u;
% Returns the prefix equivalent of the dpoly u.
if null u then 0 else dp!=replus dp!=2a u;
symbolic procedure dp!=2a u;
if null u then nil
else ((if bc_minus!? x then
list('minus,dp!=retimes(bc_2a bc_neg x . y))
else dp!=retimes(bc_2a x . y))
where x = dp_lc u, y = mo_2a dp_lmon u)
. dp!=2a cdr u;
symbolic procedure dp!=replus u;
if atom u then u else if null cdr u then car u else 'plus . u;
symbolic procedure dp!=retimes u;
% U is a list of prefix expressions the first of which is a number.
% The result is the prefix representation for their product.
if car u = 1 then if cdr u then dp!=retimes cdr u else 1
else if null cdr u then car u
else 'times . u;
% ----------- Printing routines for dpolys --------------
% ---- Authors: R. Gebauer, A. C. Hearn, H. Kredel ------
symbolic procedure dp_print u;
% Prints a distributive polynomial in infix form.
<< terpri(); dp_print1(u,nil); terpri(); terpri() >>;
symbolic procedure dp_print1(u,v);
% Prints a dpoly in infix form.
% U is a distributive form. V is a flag which is true if a term
% has preceded current form.
if null u then if null v then print_lf 0 else nil
else begin scalar bool,w;
w := dp_lc u;
if bc_minus!? w then <<bool := t; w := bc_neg w>>;
if bool then print_lf " - " else if v then print_lf " + ";
( if not bc_one!? w or mo_zero!? x then
<< bc_prin w; mo_prin(x,t)>>
else mo_prin(x,nil))
where x = dp_lmon u;
dp_print1(cdr u,t)
end;
symbolic procedure dp_print2 u;
% Prints a dpoly with restricted number of terms.
(if c and (length u>c) then
begin scalar i,v,x;
v:=for i:=1:c collect <<x:=car u; u:=cdr u; x>>;
dp_print1(v,nil); write" + # ",length u," terms #"; terpri();
end
else << dp_print1(u,nil); terpri() >>)
where c:=get('cali,'printterms);
% -------------- Auxiliary dpoly operations -------------------
symbolic procedure dp_ecart p;
% Returns the ecart of the dpoly p.
if null p then 0 else (dp!=ecart p) - (mo_ecart dp_lmon p);
symbolic procedure dp!=ecart p;
if null p then 0
else max2(mo_ecart dp_lmon p,dp!=ecart cdr p);
symbolic procedure dp_homogenize(p,x);
% Homogenize (according to mo_ecart) the dpoly p using the variable x.
if null p then p
else begin integer maxdeg;
maxdeg:=0;
for each y in p do maxdeg:=max2(maxdeg,mo_ecart car y);
return dp!=compact dp_neworder for each y in p collect
mo_inc(car y,x,maxdeg-mo_ecart car y) . cdr y;
end;
symbolic procedure dp_seed(p,s);
% Returns the dpoly p with all vars outside the list s set equal to 1.
if null p then p
else dp!=compact dp_neworder
for each x in p collect mo_seed(car x,s).cdr x;
symbolic procedure dp!=compact p;
% Collect equal terms in the sorted dpoly p.
if null p then p else dp_sum(list car p,dp!=compact cdr p);
symbolic procedure dp_xlt(p,x);
% x is the main variable. Returns the leading term of p wrt. x or p,
% if p is free of x.
if null p then p
else begin scalar d,m;
d:=mo_varexp(x,dp_lmon p);
if d=0 then return p;
return for each m in p join
if mo_varexp(x,car m)=d then {mo_inc(car m,x,-d) . cdr m};
end;
% -- dpoly operations based on computation with ideal bases.
symbolic procedure dp_pseudodivmod(g,f);
% Returns a dpoly list {q,r,z} such that z * g = q * f + r and
% z is a dpoly unit. Computes redpol({[f.e_1]},[g.0]).
% g, f and r must belong to the same free module.
begin scalar u;
f:=list bas_make1(1,f,dp_from_ei 1);
g:=bas_make(0,g);
u:=red_redpol(f,g);
return {dp_neg dp_deletecomp bas_rep car u,bas_dpoly car u,cdr u};
end;
symbolic operator dpgcd;
symbolic procedure dpgcd(u,v);
if !*mode='algebraic then dp_2a dpgcd!*(dp_from_a u,dp_from_a v)
else dpgcd!*(u,v);
symbolic procedure dpgcd!*(u,v);
% Compute the gcd of two polynomials by the syzygy method :
% 0 = u*u1 + v*v1 => gcd = u/v1 = -v/u1 .
if dp_unit!? u or dp_unit!? v then dp_fi 1
else begin scalar w;
w:=bas_dpoly first dpmat_list
syzygies!* dpmat_make(2,0,{bas_make(1,u),bas_make(2,v)},nil,nil);
return car dp_pseudodivmod(u,dp_comp(2,w));
end;
endmodule; % dpoly
end;