module kuechl;% Walking faster,B . Amrhrein,O . Gloor,W . Kuechlin
% in: Calmet,Limongelli(Eds .)Design and
% Implementation of Symbolic Computation Systems,Sept.1996
% Version 3 with a rational local solution(after letters from H.M.Moeller).
% Version 4 with keeping the polynomials as DIPs converting only
% their order mode.
switch trgroeb;
put('groebner_walk,' psopfn,' groeb!-walk);
symbolic procedure groeb!-walk u;
begin if !*groebopt then
rerror(groebner,31,"don't call 'groebner_walk' with 'on groebopt'");
if null dipvars!* then rerror(groebner,30,"'torder' must be called before");
groetime!*:=time();
!*gsugar:=t;!*groebrm:=nil;u:=car groeparams(u,1,1);
groebnervars(u,nil);u:=groeb!-list(u,'simp);
groedomainmode();u:=groeb!-w2 u;
return'list.groeb!-collect(u,'mk!*sq)end;
symbolic procedure groeb!-list(u,fcn);
% Execute the function ' fcn ' for the elements of the algebriac
% list 'u'.
<<if atom u or not(eqcar(u,'list))then
rerror('groebner,29,"groebner: list as argument required");
groeb!-collect(cdr u, fcn)>>;
symbolic procedure groeb!-collect(l,f);
% Collect the elements of function 'f' applied to the elements of
% the symbolic list 'l'. If 'f' is a number,map 'l' only.
for each x in l collect if numberp f then f else apply1(f,x);
symbolic procedure groeb!-w2 g;
% This is(essentially)the routine Groebner_Walk.
% G is a list of standard quotients,
% a Groebner basis gradlex or based on a vector like [1 1 1 ...].
% The result is the Groebner basis(standard quotients)with the
% final term order(lex)as its main order.
begin scalar iwv,owv,omega,gomega,gomegaplus,tt,tto,pc;
scalar first,mx,imx,mmx,immx,nn,ll,prim;
scalar !*vdpinteger,!*groebdivide;
!*vdpinteger: nil; % switch on division mode
!*groebdivide:=t;
first:=t;pcount!*:=0;mmx:=!*i2rn 1;immx:=mmx;
iwv:=groeb!-collect(dipvars!*,1);
omega:=iwv; % Input order vector.
owv:=1 .groeb!-collect(cdr dipvars!*,0);
tto:=owv; % Output order vector .
groeb!-w9('weighted,omega);% Install omega as weighted order.
g:=groeb!-collect(g,'sq2vdp);
pc:=pcount!*;
gbtest g; % Test the Groebner property.
nn:=length dipvars!*;
ll:=rninv!: !*i2rn nn; % Inverse of the length.
prim:=t; % Preset.
loop:groeb!-w9(' weighted,omega);
mx:=groeb!-w6!-4 groeb!-collect(omega,1);
% Compute the maximum of \omega.
if !*trgroeb then groebmess34 cadr mx;
imx:=rninv!: mx;
g:=if first then groeb!-collect(g,'vdpsimpcont) else groeb!-w10 g;
if !*trgroeb then groebmess29 omega;
gomega:=if first or not prim then g
else groeb!-w3(g,omega); % G_\omega = initials(G_\omega);
pcount!*:=pc;
if !*trgroeb and not first then groebmess32 gomega;
gomegaplus:=if first then{gomega}else gtraverso(gomega,nil,nil);
if cdr gomegaplus then rerror(groebner,31,
"groebner_walk,cdr of 'groebner' must be nil")
else gomegaplus:=car gomegaplus;
if !*trgroeb and not first then groebmess30 gomegaplus;
if not first and prim
then g:=groeb!-w4(gomegaplus,gomega,g)
else if not prim then g:=gomega;
% G=lift(G_{%omega}{plus},<{plus},G_{%omega),G,<)
if not first then g:=for each x in g collect gsetsugar(x,nil);
if !*trgroeb and not first then groebmess31 g;
if groeb!-w5(omega,imx,tto,immx)then go to ret;
% Stop if tt has been 1 once.
if not first and rnonep!: tt then go to ret;% Secodary abort crit.
tt:=groeb!-w6!-6(g,tto,immx,omega,imx,ll);% Determine_border .
if !*trgroeb then groebmess36 tt;
if null tt then go to ret;
% criterion: take primary only if tt neq 1
prim:=not rnonep!: tt;
if !*trgroeb then groebmess37 prim;
%\omega =(1-t)*\omega+t*tau
omega:=groeb!-w7(tt,omega,imx,tto,immx);
if !*trgroeb then groebmess35 omega;
first:=nil;go to loop;
ret:
if !*trgroeb then groebmess33 g;
g:=groeb!-collect(g,'vdpsimpcont);
g:=groeb!-collect(g,'vdp2sq);
return g end;
symbolic procedure groeb!-w3(g,omega);
% Extract head terms of g corresponding to omega.
begin scalar x,y,gg,ff;
gg:=for each f in g collect<<ff:=vdpfmon(vdplbc f,vdpevlmon f);
gsetsugar(ff,nil);
x:=evweightedcomp2(0,vdpevlmon ff,omega);
y:=x;
f:=vdpred f;
while not vdpzero!? f and y=x do
<<y:=evweightedcomp2(0,vdpevlmon f,omega);
if y=x then
ff:=vdpsum(ff,vdpfmon(vdplbc f,vdpevlmon f));f:=vdpred f>>;
ff>>;return gg end;
symbolic procedure groeb!-w4(gb,gomega,g);
% gb Groebner basis of gomega,
% gomega head term system g_\omega of g,
% g full(original)system of polynomials.
begin scalar x;
for each y in gb do gsetsugar(y,nil);
x:=for each y in gomega collect groeb!-w8(y,gb);
x:=for each z in x collect groeb!-w4!-1(z,g);return x end;
symbolic procedure groeb!-w4!-1(pl,fs);
% pl is a list of polynomials corresponding to the full system fs.
% Compute the sum of pl*fs. Result is the sum.
begin scalar z;z:=vdpzero();
gsetsugar(z,0);
for each p in pair(pl,fs)do
if car p then z:=vdpsum(z,vdpprod(car p,cdr p));
z:=vdpsimpcont z;return z end;
symbolic procedure groeb!-w5(ev1,x1,ev2,x2);
% ev1=ev2 equality test.
groeb!-w5!-1(x1,ev1,x2,ev2);
symbolic procedure groeb!-w5!-1(x1,ev1,x2,ev2);
( null ev1 and null ev2)or
(rntimes!:(!*i2rn car ev1,x1)=rntimes!:(!*i2rn car ev2,x2)
and groeb!-w5!-1(x1,cdr ev1,x2,cdr ev2));
symbolic procedure groeb!-w6!-4 omega;
% Compute the weighted length of \omega.
groeb!-w6!-5(omega,vdpsortextension!*,0);
symbolic procedure groeb!-w6!-5(omega,v,m);
if null omega then !*i2rn m
else if 0=car omega then groeb!-w6!-5(cdr omega,cdr v,m)
else if 1 = car omega
then groeb!-w6!-5(cdr omega,cdr v,m #+ car v)
else groeb!-w6!-5(cdr omega,cdr v,m #+ car omega #* car v);
symbolic procedure groeb!-w6!-6(gb,tt,ifactt,tp,ifactp,ll);
% Compute the weight border(minimum over all polynomials of gb).
begin scalar mn,x,zero,one;
zero:=!*i2rn 0;one:=!*i2rn 1;
while not null gb do
<<x:=groeb!-w6!-7(car gb,tt,ifactt,tp,ifactp,zero,one,ll);
if null mn or(x and rnminusp!: rndifference!:(x,mn))
then mn:=x;gb:=cdr gb>>;return mn end;
symbolic procedure groeb!-w6!-7(pol,tt,ifactt,tp,ifactp,zero,one,ll);
% Compute the minimal weight for one polynomial;the idea is,
% that the polynomial has a degree greater than 0.
begin scalar a,b,ev1,ev2,x,y,z,mn;
ev1:=vdpevlmon pol;
a:=evweightedcomp2(0,ev1,vdpsortextension!*);
y:=groeb!-w6!-8(ev1,tt,ifactt,tp,ifactp,zero,zero,one,ll);
y:=(rnminus!: car y).(rnminus!: cdr y);
pol:=vdpred pol;
while not(vdpzero!? pol)do
<<ev2:=vdpevlmon pol;
pol:=vdpred pol;
b:=evweightedcomp2(0,ev2,vdpsortextension!*);
if not(a=b)then
<<x:=groeb!-w6!-9(ev2,tt,ifactt,tp,ifactp,car y,cdr y,one,ll,nil);
if x then
<<z:=rndifference!:(x,one);
if rnminusp!: rndifference!:(zero,x)and
(rnminusp!: z or rnzerop!: z)and
(null mn or rnminusp!: rndifference!:(x,mn))
then mn:=x>>>>>>;return mn end;
symbolic procedure groeb!-w6!-8(ev,tt,ifactt,tp,ifactp,sum1,sum2,m,dm);
begin scalar x,y,z;
if ev then<<x:=rntimes!:(!*i2rn car ev,m);
y:=rntimes!:(!*i2rn car tp,ifactp);
z:=rntimes!:(!*i2rn car tt,ifactt)>>;
return
if null ev then sum1.sum2 else
groeb!-w6!-8(cdr ev,cdr tt,ifactt,cdr tp,ifactp,
rnplus!:(sum1,rntimes!:(y,x)),
rnplus!:(sum2,rntimes!:(rndifference!:(z,y),x)),
rndifference!:(m,dm),
dm)end;
symbolic procedure groeb!-w6!-9(ev,tt,ifactt,tp,ifactp,y1,y2,m,dm,done);
% Compute the rational solution s:
%(tp+s*(tt-tp))*ev1=(tp+s*(tt-tp))*evn.
% The sum with ev1 is collected already in y1 and y2(with negative sign).
% This routine collects the sum with evn and computes the solution.
begin scalar x,y,z;
if ev then<<x:=rntimes!:(!*i2rn car ev,m);
y:=rntimes!:(!*i2rn car tp,ifactp),
z:=rntimes!:(!*i2rn car tt,ifactt)>>;
return if null ev then
if null done then nil
else rnquotient!:(rnminus!: y1,y2)
else
groeb!-w6!-9(cdr ev,cdr tt,ifactt,cdr tp,ifactp,
rnplus!:(y1,rntimes!:(y,x)),
rnplus!:(y2,rntimes!:(rndifference!:(z,y),x)),
rndifference!:(m,dm),
dm,
done or not(car ev = 0)) end;
symbolic procedure groeb!-w7(tt,omega,x,tto,y);
% Compute omega*x*(1-tt)+tto*y*tt.
% tt is a rational number.
% x and y are rational numbers(inverses of the legths of omega/tt).
begin scalar n,z;n:=!*i2rn 1;
omega:=for each g in omega collect
<<z:=rnplus!:(rntimes!:(rntimes!:(!*i2rn g,x),
rndifference!:(!*i2rn 1,tt)),
rntimes!:(rntimes!:(!*i2rn car tto,y),tt));
tto:=cdr tto;
n:=groeb!-w7!-1(n,rninv!: z);z>>;
omega:=for each a in omega collect rnequiv rntimes!:(a,!*i2rn n);
return omega end;
symbolic procedure groeb!-w7!-1(n,m);
% Compute lcm of n and m. N and m are rational numbers.
% Return the lcm.
% Ignore the denominators of n and m.
begin scalar x,y,z;
if atom n then x:=n else
<<x:=rnprep!: n;
if not atom x then x:=cadr x>>;
if atom m then y:=m else
<<y:=rnprep!: m;if not atom y then y:=cadr y>>;
z:=lcm(x,y);return z end;
symbolic procedure groeb!-w8(p,gb);
% Computes the cofactor of p wrt gb.
% Result is a list of cofactors corresponding to g.
% The cofactor 0 is represented as nil.
begin scalar x,y;
x:=groeb!-w8!-1(p,gb);p:=secondvalue!*;
while not vdpzero!? p do
<<y:=groeb!-w8!-1(p,gb);p:=secondvalue!*;
x:=for each pp in pair(x,y)do
if null car pp then cdr pp
else if null cdr pp then car pp
else vdpsum(car pp,cdr pp)>>;return x end;
symbolic procedure groeb!-w8!-1(p,gb);
% Search in groebner basis gb the polynomial which divides the
% head monomial of the polynomial p. The walk version of
% groebsearchinlist.
% Result: the sequence corresponding to g with the monomial
% factor inserted.
begin scalar e,cc,r,done,pp;
pp:=vdpevlmon p;cc:=vdplbc p;
r:=for each poly in gb collect
if done then nil else
if vevdivides!?(vdpevlmon poly,pp)then
<<done:=t;e:=poly;
cc:=vbcquot(cc,vdplbc poly);
pp:=vevdif(pp,vdpevlmon poly);
secondvalue!*:=
vdpsum(vdpprod(gsetsugar(vdpfmon(vbcneg cc,pp),nil),poly), p);
vdpfmon(cc,pp)>>
else nil;
if null e then
<<print p;print "-----------------";print gb;
rerror(groebner,28,"groeb-w8-1 illegal structure")>>;
return r end;
symbolic procedure groeb!-w9(mode,ext);
% Switch on vdp order mode 'mode' with extension 'ext'.
% Result is the previous extension.
begin scalar x;
x:=vdpsortextension!*;vdpsortextension!*:=ext;
torder2 mode;return x end;
symbolic procedure groeb!-w10 s;
% Convert the dips in s corresponding to the actual order.
groeb!-collect(s,'groeb!-w10!-1);
symbolic procedure groeb!-w10!-1 p;
% Convert the dip p corresponding to the actal order.
begin scalar x;
x:=vdpfmon(vdplbc p,vdpevlmon p);
x:=gsetsugar(vdpenumerate x,nil);
p:=vdpred p;
while not vdpzero!? p do
<<x:=vdpsum(vdpfmon(vdplbc p,vdpevlmon p),x);p:=vdpred p>>;
return x end;
symbolic procedure rninv!: x;
% Return inverse of a(rational)number x: 1/x.
<<if atom x then x:=!*i2rn x;
x:=cdr x;
if car x < 0 then mkrn(- cdr x,- car x)else mkrn(cdr x,car x)>>;
symbolic procedure sq2vdp s;
% Standard quotient to vdp.
begin scalar x,y;x:=f2vdp numr s;
gsetsugar(x,nil);y:=f2vdp denr s;
gsetsugar(y,0);s:=vdpdivmon(x,vdplbc y,vdpevlmon y);
return s end;
symbolic procedure vdp2sq v;
% Conversion vdp to standard quotient.
begin scalar x,y,z,one;one := 1 ./ 1;x := nil ./ 1;
while not vdpzero!? v do
<<y:=vdp2f vdpfmon(one,vdpevlmon v)./ 1;
z:=vdplbc v;
x:=addsq(x,multsq(z,y));
v:=vdpred v>>;return x end;
endmodule;;end;