module invbase; % Computing involutive basis of polynomial system.
% Authors: Alexey Yu. Zharkov, Yuri A. Blinkov
% Saratov University, Astrakhanskaya 83,
% Saratov 410071, Russia
% e-mail: postmaster@scnit.saratov.su
% Copyright A. Zharkov, Y. Blinkov;
% all rights reserved.
% Minor fixes by John Fitch.
create!-package('(invbase invbint invbcomp),'(contrib invbase));
fluid '(conds gv hv bv ng gg varlist vjets nc); % globals
fluid '(ordering redtails); % modes
fluid '(path tred stars); % tracing
fluid '(reductions nforms zeros maxord title); % statistics
fluid '(invsysvars!* !*trinvbase alfa beta shortway thirdway
invtempbasis);
share invtempbasis;
ordering := 'grev;
switch trinvbase;
gv:=mkvect(1000)$ % p o l y n o m i a l s
bv:=mkvect(1000)$ % f l a g (n e w p r o l o n g a t i o n s)
endmodule;
module invbint; % Algebraic mode interface to invbase.
symbolic procedure invtorder u;
begin scalar w,o;
w := reval car u;
o := assoc(w,'((gradlex . glex) (revgradlex .grev) (lex . lex)));
if null o then typerr(w,"involutive term ordering");
ordering := cdr o;
invsysvars!* :=
if cdr u then
for each y in cdr listeval (cadr u,nil) collect reval y
else nil;
end;
put('invtorder,'stat,'rlis);
symbolic procedure invbase u;
begin scalar sys,vars,r;
u := reval car u;
if not eqcar(u,'list) then rederr "Argument to invbase not a list";
sys := for each p in cdr u collect
<<p := reval p;
if eqcar(p,'equal) then p:=reval{'difference,cadr p,caddr p};
p>>;
% find the variables.
vars := invsysvars!* or gvarlis sys;
readsys('list.sys,'list.vars);
invbase!*();
r:= for each p in gg collect
'plus . for each m in getv(gv,car p) collect
prepsq !*di2q(list m,vars);
return 'list . r;
end;
put('invbase,'psopfn,'invbase);
symbolic procedure invlex u;
begin scalar sys,vars,r;
u := reval car u;
if not eqcar(u,'list) then rederr "Argument to invlex not a list";
sys := for each p in cdr u collect
<<p := reval p;
if eqcar(p,'equal) then p:=reval{'difference,cadr p,caddr p};
p>>;
% find the variables.
vars := invsysvars!* or gvarlis sys;
readsys('list.sys,'list.vars);
invlex!*();
(r:= for each p in gg collect
'plus . for each m in getv(gv,car p) collect
prepsq !*di2q(list m,vars)) where ordering='lex;
return 'list . r;
end;
put('invlex,'psopfn,'invlex);
symbolic procedure invtest u;
begin scalar sys,vars,r;
u := reval car u;
if not eqcar(u, 'list) then rederr "Argument to invtest not a list";
sys := for each p in cdr u collect
<<p := reval p;
if eqcar(p,'equal) then p:=reval{'difference,cadr p,caddr p};
p>>;
% find the variables.
vars := invsysvars!* or gvarlis sys;
readsys('list.sys,'list.vars);
return invtest!*();
end;
put('invtest,'psopfn,'invtest);
% the following procedure are borrowed from the groebner package:
symbolic procedure gvarlis u;
% Finds variables (kernels) in the list of expressions u.
sort(gvarlis1(u,nil),function ordop);
symbolic procedure gvarlis1(u,v);
if null u then v
else union(gvar1(car u,v),gvarlis1(cdr u,v));
symbolic procedure gvar1(u,v);
if null u or numberp u or (u eq 'i and !*complex) then v
else if atom u then if u member v then v else u . v
else if get(car u,'dname) then v
else if car u memq '(plus times expt difference minus)
then gvarlis1(cdr u,v)
else if car u eq 'quotient then gvar1(cadr u,v)
else if u member v then v
else u . v;
endmodule;
module invbcomp;
%----------------------------------------------------------------------
symbolic procedure c_zero(); nil$ % REPRESENTATION OF ZERO
%----------------------------------------------------------------------
symbolic procedure cneg(c); % - C
negf c$
%----------------------------------------------------------------------
symbolic procedure csum(c1,c2); % C1 + C2
addf(c1,c2);
%----------------------------------------------------------------------
symbolic procedure cprod(c1,c2); % C1 * C2
multf(c1,c2);
%----------------------------------------------------------------------
symbolic procedure cdiv(c1,c2); % C1/C2
numr resimp(c1 . c2);
%----------------------------------------------------------------------
symbolic procedure trass(id,value); % tracing of assignments
<< terpri(); write id; write " = "; write value; terpri(); >>$
%----------------------------------------------------------------------
symbolic procedure leftzeros(u); % u : list
if null u or car u neq 0 then 0 else 1 #+ leftzeros cdr u$
%----------------------------------------------------------------------
procedure class(jet);
if ord jet = 0 then 0 else 1
#+ leftzeros reverse (if ordering = 'lex then jet else cdr jet);
%----------------------------------------------------------------------
symbolic procedure ord(jet);
if ordering = 'lex then eval('plus . jet) else car jet$
%----------------------------------------------------------------------
symbolic procedure ljet(p); caar p$
%----------------------------------------------------------------------
symbolic procedure sub01(v,u);
%%% replace each x in u by < if x=v then 1 else 0 >
if u then (if (car u = v) then 1 else 0) . sub01(v,cdr u);
%----------------------------------------------------------------------
symbolic procedure !*v2j(v);
if ordering = 'lex then sub01(v,varlist) else (1 . sub01(v,varlist) );
%----------------------------------------------------------------------
symbolic procedure nonmult(cl); % --> list of vjets
reverse cdr member(nth(reverse vjets,cl),reverse vjets);
%----------------------------------------------------------------------
symbolic procedure insert(x,gg);
begin scalar gg1;
while gg and dless(cdr x,cdar gg) do
<< gg1 := car gg . gg1; gg := cdr gg >>;
return append(reversip gg1, x . gg);
end;
%----------------------------------------------------------------------
symbolic procedure addnew(f,ind,ff);
%%% adds element f to set (with index ind), returns modified ff
<< putv(gv,ind,f); putv(bv,ind,t);
if null f then ff
else ff := insert(ind . ljet f, ff)
>>$
%----------------------------------------------------------------------
symbolic procedure dlesslex(d1,d2);
%%%% RETURNS T IF D1 < D2 (lex), NIL OTHERWISE
if null d1 then nil else
if car d1 #> car d2 then nil else
if car d1 #< car d2 then t else dlesslex(cdr d1,cdr d2);
%----------------------------------------------------------------------
symbolic procedure dless(d1,d2); % --> T if d1 < d2 , NIL otherwise
if ordering = 'lex then dlesslex(d1,d2) else
if car d1 #< car d2 then t else if car d1 #> car d2 then nil
else if ordering = 'glex then dlesslex(cdr d1,cdr d2)
else if ordering = 'grev then dlesslex(reverse cdr d2, reverse cdr d1);
%-----------------------------------------------------------------------
symbolic procedure ddmult(d1,d2);
if null d1 then nil else (car d1 #+ car d2) . ddmult(cdr d1,cdr d2);
%-----------------------------------------------------------------------
symbolic procedure dquot(d2,d1);
%%%% RETURNS D2-D1 IF D1 DIVIDES D2, NIL OTHERWISE
begin scalar d3; integer n;
l1:n:=car(d2)-car(d1);
if n #< 0 then return nil;
d3:=n . d3;
d1:=cdr d1; d2:=cdr d2;
if d1 then goto l1;
return reversip d3;
end;
%-----------------------------------------------------------------------
symbolic procedure pcmult(p,c); % P*C (C IS NOT ZERO)
for each x in p collect car x.cprod(c,cdr x);
%-----------------------------------------------------------------------
symbolic procedure pcdiv(p,c); % P/C (division in ring)
for each x in p collect car x . cdiv(cdr x,c);
%-----------------------------------------------------------------------
symbolic procedure pdmult(p,d); % P*< D >
for each x in p collect
(for each y in pair(car x,d) collect car(y)#+cdr(y)).cdr x$
%-----------------------------------------------------------------------
symbolic procedure psum(p1,p2); % P1 + P2
begin scalar t1,t2,d2,c3,p3,sum,ret;
if null p1 then sum:=p2 else
if null p2 then sum:=p1 else
while p2 and not ret do
<< t2:=car p2; d2:=car t2;
while p1 and dless(d2,caar p1) do
<< p3:=car(p1).p3; p1:=cdr p1 >>;
if null p1 then
<< sum:=append(reverse p3,p2); ret:=t >> else
<< t1:=car p1;
if d2=car t1 then %%%% NOW T1<=T2
<< c3:=csum(cdr t1,cdr t2); %%%% LIKE TERM
if c3 neq c_zero() then p3:=(d2.c3).p3;
p1:=cdr p1;
t1:=if p1 then car p1; %%%% NEW T1
>>
else p3:=t2.p3; %%%% OLD T1
p2:=cdr p2; %%%% NEW T2
if null p2 then sum:= append(reverse p3,p1)
>> >>;
return sum
end;
%-----------------------------------------------------------------------
symbolic procedure pneg(p); % - P
for each x in p collect car(x).cneg(cdr(x));
%-----------------------------------------------------------------------
symbolic procedure pdif(p1,p2); % P1 - P2
psum(p1,pneg p2);
%-----------------------------------------------------------------------
symbolic procedure dd(d1,d2); % uses fluid VJETS
begin scalar dq,lz;
dq:=dquot(d2,d1);
if not dq then return if dless(d1,d2) then 1 % D1 < D2
else 0; % D1 > D2
if ordering neq 'lex then dq:=cdr dq;
lz := leftzeros dq;
return
if not nc and not(lz #< length varlist #- class d1)
then 3 % D1 divides D2 (mult.)
else if nc and not(lz #< length varlist #- nc)
then 4 % D1 divides D2 in 1:nc classes and coincides in others
else 2; % D1 divides D2 (usual)
end;
%-----------------------------------------------------------------------
symbolic procedure dlcm(d1,d2);
if ordering='lex then for each x in pair(d1,d2) collect max(car x,cdr x)
else addgt( for each x in pair(cdr d1,cdr d2) collect max(car x,cdr x));
%-----------------------------------------------------------------------
symbolic procedure nf(h,gg,sw);
%%%% H = NORMALIZED POLYNOMIAL
%%%% GG = LIST OF KEYED LPP'S OF GG-SET
%%%% RETURNS NORMAL FORM OF H WITH RESPECT TO GG-SET
%%%% ===============================================
if null gg then h else
begin scalar f,lpf,g,c,cf,cg,nf,g1,g2,u,nr; nr:=0;
f:=h; g1:=gg;
nextlpf: if null f then goto exit;
lpf:=caar f;
% diminish G1 so that LPF >= G1 (and might be reduced !)
% ------------------------------------------------------
while not null g1 and dless(lpf,cdar g1) do g1:=cdr g1;
if null g1 then goto exit;
g2:=g1; % NOW LPF >= G2
% reduction of LPF
% ----------------
while g2 and dd(cdar g2,lpf) #< sw + 2 do g2:=cdr g2;
if null g2 then % LPF NOT REDUCED
( if redtails then << nf:=(lpf.cdar f).nf; f:=cdr f >>
else goto exit )
else % REDUCTION OF LPF
<< g:=getv(gv,caar g2);
c:=gcdf!*(cdar f, cdar g);
cf:=cdiv(cdar f,c); cg:=cdiv(cdar g,c);
f:=pcmult(f,cg); nf:=pcmult(nf,cg); g:=pcmult(g,cf);
u:=pdmult(cdr g, dquot(lpf,cdar g2));
if tred then
<< terpri();
write "r e d u c t i o n : ",lpf,"/",cdar g2;
terpri();
>>;
if stars then write "*";
nr := nr #+ 1;
f:=pdif(cdr f,u);
>>;
goto nextlpf;
exit:
reductions := reductions #+ nr;
nforms := nforms #+ 1;
u:= gcdout append(reversip nf,f);
if null u then zeros := zeros #+ 1;
return u;
end;
%-----------------------------------------------------------------------
symbolic procedure gcdout(p);
% cancel coeffs of P by their common factor.
if !*modular then p else
if null p then nil else if ord ljet p = 0 then p else
begin scalar c,p1;
p1:=p; c:=cdar p1;
while p1 and c neq 1 do << c:=gcdf!*(c,cdar p1); p1:=cdr p1 >>;
return if c = 1 then p else pcdiv(p,c);
end;
%-----------------------------------------------------------------------
expr procedure newbasis(gg,sw)$
%%%% SIDE EFFECT: CHANGES CDR'S OF GV(K);
begin scalar g1,g2;
g1:=reverse gg;
while g1 do
<< putv(gv,caar g1,nf(getv(gv,caar g1),g2,sw));
g2:=(car g1).g2; g1:=cdr g1;
>>;
end$
%-----------------------------------------------------------------------
symbolic procedure !*f2di(f,varlist);
%%% f: st.f., varlist: kernel list --> f in distributive form
if null f then nil else
if domainp f then
((addgt for each v in varlist collect 0).(f)).nil else
psum(if member(mvar f,varlist) then
pdmult(!*f2di(lc f,varlist),
addgt for each v in varlist collect
if v = mvar f then ldeg f else 0
)
else pcmult(!*f2di(lc f,varlist),((lpow f.1).nil)),
!*f2di(red f,varlist) );
%-----------------------------------------------------------------------
symbolic procedure !*di2q0(p,varlist);
if null p then nil . 1 else
addsq( (lambda s,u;
<< for each x in u do
if cdr x neq 0 then s:=multsq(s,((x.1).nil).1);
s >>
) (cdar p, pair(varlist,
if ordering='lex then ljet p else cdr ljet p)),
!*di2q0(cdr p,varlist) );
%----------------------------------------------------------------------
symbolic procedure !*di2q(p,varlist);
!*di2q0(for each x in p collect car x . (cdr x . 1), varlist);
%----------------------------------------------------------------------
symbolic procedure show(str,p); % p = poly in a special (dist.) form
if null p then (algebraic write str," := 0")
else algebraic write str," := ",
lisp prepsq !*di2q(list car p, varlist)," + ",
lisp prepsq !*di2q(cdr p, varlist);
%----------------------------------------------------------------------
lisp procedure addgt(u);
if ordering = 'lex then u else eval('plus.u) . u$
%-----------------------------------------------------------------------
symbolic procedure printsys(str,gg);
begin scalar i; i:=0;
for each x in gg do
<< i:=i+1;
algebraic write str,"(",lisp i,") := ",
lisp prepsq !*di2q(list car getv(gv,car x), varlist)," + ",
lisp prepsq !*di2q(cdr getv(gv,car x), varlist);
>>;
end;
%-----------------------------------------------------------------------
symbolic procedure answer(gg);
<< if title then algebraic write "% ",lisp title;
trass("% ORDERING",varlist); printsys("G",reverse gg);
>>$
%-----------------------------------------------------------------------
symbolic procedure wr(file,gg);
<< off nat,time; out file;
write "algebraic$"; write "operator g$";
answer(gg);
write "end;"; shut file; on nat,time >>$
%-----------------------------------------------------------------------
symbolic procedure invtest!*();
begin scalar g,c; c:=t;
if !*trinvbase then terpri();
for each x in gg do
if c then
<< g:=getv(gv, car x);
for each vj in nonmult(class ljet g) do
if c and nf(pdmult(g,vj),gg,1) then
<< c:=nil;
if !*trinvbase then prin2t "INV - t e s t f a i l e d"; >>;
>>;
if c and !*trinvbase then prin2t "I n v o l u t i v e b a s i s";
return c;
end;
%-----------------------------------------------------------------------
symbolic procedure redall(gg,ff,sw);
% side effect : changes flag thirdway.
begin scalar rr,f,f1,lj,k,new;
rr := ff; thirdway:=shortway:=nil; new:=t;
while rr do
<< f:=car reverse rr; rr:=delete(f,rr);
k:=car f; f1:=getv(gv,k);
if path then
<< % write k,": ";
if new then write ljet f1," ==> "
else write ljet f1," --> ";
>>;
f:=putv(gv,k,nf(f1,gg,sw));
lj:=if f then ljet f else 0;
if path then
<< write lj; terpri() >>;
if null f then nil else
if ord lj = 0 then conds := f . conds else
<< if ljet f neq ljet f1 then shortway:=t;
if not new and f neq f1 then thirdway:=t;
for each x in gg do if dd(lj,cdr x) >= sw + 2 then
<< gg:=delete(x,gg); rr:=insert(x,rr);
putv(bv,car x,t); %
>>;
gg:=insert(k.lj,gg); new:=nil;
>> >>;
return gg;
end;
%-----------------------------------------------------------------------
symbolic procedure remred(ff,sw); % removes redundant elements
begin scalar gg,gg1,f,g,p;
ff:=reverse ff;
while ff do
<< f:=car ff; ff:=cdr ff;
p:=t; gg1:=gg;
while p and gg1 do
<< g:=car gg1; gg1:=cdr gg1;
if dd(cdr g,cdr f) >= sw + 2 then p:=nil;
>>;
if p then gg := f . gg;
>>;
return gg;
end;
%-----------------------------------------------------------------------
symbolic procedure invbase!*();
begin scalar gg1,g,k,nm,f,thirdway,shortway,fin,p,p0,lb,r;
if !*trinvbase then terpri();
p:=maxord:=-1;
if path then terpri();
gg:=redall(nil,gg,1);
newbasis(gg,1);
lb:=0;
for each x in gg do lb:=lb + ord cdr x;
lb:=lb + length varlist - 1;
l: gg1 := reverse gg;
while gg1 and null getv(bv,caar gg1) do gg1 := cdr gg1;
if gg1 then
<< if cadar gg1 = cadar gg then
<< p0:=p;
p:=cadar gg1;
if !*trinvbase and p > p0 then
<< terpri();
write "---------- ORDER = ",cadar gg," ----------";
terpri(); terpri();
>>;
if p > lb then
<< gg:=redall(nil,gg,0);
newbasis(gg,0);
invtempbasis := 'list .
for each x in gg
collect 'plus . for each m in getv(gv,car x)
collect prepsq !*di2q(list m,varlist);
rederr "Maximum degree bound exceeded.";
>>;
maxord:=max(maxord,cadar gg);
if cadar gg < maxord then fin:=t;
>>;
if fin then goto m;
k := caar gg1;
g := getv(gv,k); putv(bv,k,nil);
nm := nonmult(class ljet g);
for each vj in nm do
<< ng := ng + 1;
f := pdmult(g,vj); putv(gv,ng,f); putv(bv,ng,t);
gg := redall(gg,list(ng.ljet f),1);
if thirdway then newbasis(gg,1) else
if shortway then for each y in gg do if car y neq ng then
putv(gv,car y,nf(getv(gv,car y),list(ng.ljet getv(gv,ng)),1));
>>;
go to l;
>>;
m: stat(); if p <= lb then dim gg;
end;
%-----------------------------------------------------------------------
symbolic procedure njets(n,q); % number of jets of n vars and order q
combin(q,q+n-1);
%----------------------------------------------------------------------
symbolic procedure combin(m,n); % number of combinations of m from n
if m>n then 0 else
begin integer i1,i2; i1:=i2:=1;
for i:=n-m+1:n do i1:=i1*i; for i:=1:m do i2:=i2*i;
return i1/i2;
end;
%----------------------------------------------------------------------
symbolic procedure dim(gg);
if !*trinvbase then
begin integer q,n,cl,s,y,dim,dp,mon;
q:=cadar gg; n:=length varlist; dim:=0;
for i:=1:n do putv(beta,i,0);
for each x in gg do
<< cl:=class cdr x;
for i:=cl step -1 until 1 do
<< y:=njets(cl-i+1,q-ord cdr x);
putv(beta,i,getv(beta,i)+y);
>> >>;
terpri();
for i:=1:n do
<< putv(alfa,i,combin(n-i,q+n-i)-getv(beta,i));
if getv(alfa,i) neq 0 then dim := dim + 1;
% write "a[",i,",",q,"]=",getv(alfa,i)," ";
>>;
terpri();
terpri(); write "D i m e n s i o n = ",dim; terpri();
if dim = 0 then nroots gg;
end;
%----------------------------------------------------------------------
symbolic procedure nroots(gg);
% number of roots of zero-dimensional Ideal.
if gg then begin integer d;
for each x in gg do d := d + car reverse x;
terpri(); write "N u m b e r o f s o l u t i o n s = ",d;
terpri();
end;
%----------------------------------------------------------------------
symbolic procedure stat();
if !*trinvbase then
<< terpri();
write "reductions = ",reductions;
write " zeros = ",zeros; write " maxord = ",maxord;
write " order = ",cadar gg; write " length = ",length gg;
>>$
%----------------------------------------------------------------------
symbolic procedure !*g2lex(p);
% works correctly only when ORDERING= lex.
%%% p: poly in graduate form ---> lexicographic form
begin scalar p1;
for each x in p do p1:=psum(p1,list(cdar x . cdr x));
return p1;
end;
%----------------------------------------------------------------------
symbolic procedure lexorder(lst);
% works correctly only when ORDERING = lex.
begin scalar lst1,lj;
for each x in lst do
<< lj:=ljet putv(gv, car x, gcdout !*g2lex getv(gv,car x));
lst1 := insert((car x).lj, lst1);
>>;
return lst1;
end;
%----------------------------------------------------------------------
symbolic procedure gi(gg,i); % subsystem of GG of class = i
begin scalar ff;
for each x in gg do if class cdr x = i then ff := x . ff;
return ff;
end;
%----------------------------------------------------------------------
symbolic procedure monic(jet,cl); % for lexicoraphy only
begin scalar u,n;
jet:=reverse jet;
n:=length varlist;
for i:=1:n do if i neq cl then u:=nth(jet,i).u;
return u = for each v in cdr varlist collect 0$
end;
%----------------------------------------------------------------------
symbolic procedure monicmember(gg,cl);
begin scalar p;
l: if null gg then return nil;
if monic(cdar gg,cl) then return t;
gg:=cdr gg;
go to l;
end;
%----------------------------------------------------------------------
symbolic procedure montest(gg);
begin scalar p,n;
p:=t;
n:=length varlist;
for i:=1:n do if not monicmember(gg,i) then << p:=nil; i:=n+1 >>;
return p;
end;
%----------------------------------------------------------------------
symbolic procedure invlex!*(); % side effect: changes GG
begin scalar gi,n,gginv,ordering;
n:=length varlist; gginv:=mkvect n;
ordering:='lex;
for i:=1:n do putv(gginv,i,lexorder gi(gg,i));
gg:=nil;
for i:=1:n do
<< nc:=i;
if path then << trass("i",i); terpri() >>;
gg:=redall(gg,getv(gginv,i),2);
if montest gg then i:=n+1;
>>;
nc:=nil;
gg:=remred(gg,0);
newbasis(gg,0);
end;
%----------------------------------------------------------------------
symbolic procedure readsys(elist,vlist);
begin;
varlist:=cdr vlist;
ng:=reductions:=nforms:=zeros:=0;
alfa:=mkvect length varlist;
beta:=mkvect length varlist;
gg:=nil;
for each x in cdr elist do
gg:=addnew(gcdout !*f2di(numr simp x, varlist), ng:=ng+1, gg);
vjets:=for each v in varlist collect !*v2j(v);
end;
%-----------------------------------------------------------------------
lisp operator readsys$
%-----------------------------------------------------------------------
% D E F A U L T V A L U E S
% ======================================================================
ordering:='grev$ redtails:=t$
tred := path := stars := nil$
% ======================================================================
endmodule;
end;