module glexconv; % Newbase-Algorithm:
% Faugere, Gianni, Lazard, Mora
fluid '(!*factor !*complex !*exp !*gcd !*ezgcd); % REDUCE switch
fluid '( % switches from the user interface
!*groebopt !*groebfac !*groebres !*trgroeb !*trgroebs !*groebrm
!*trgroeb1 !*trgroebr !*groebprereduce groebrestriction!*
!*fullreduction !*groebstat !*groebprot !*gltbasis
!*groebsubs
!*vdpinteger !*vdpmodular % indicating type of algorithm
!*gsugar % sugar mode
vdpsortmode!* % term ordering mode
secondvalue!* thirdvalue!* % auxiliary: multiple return values
fourthvalue!*
factortime!* % computing time spent in factoring
factorlvevel!* % bookkeeping of factor tree
pairsdone!* % list of pairs already calculated
probcount!* % counting subproblems
vbccurrentmode!* % current domain for base coeffs.
vbcmodule!* % for modular calculation: current prime
glexdomain!* % information wrt current domain
!*gsugar
);
global '(groebrestriction % interface: name of function
groebresmax % maximum number of internal results
gvarslast % output: variable list
groebprotfile
gltb
);
flag ('(groebrestriction groebresmax gvarslast groebprotfile
gltb),'share);
switch groebopt,groebfac,groebres,trgroeb,trgroebs,trgroeb1,
trgroebr,groebstat,gltbasis;
% variables for counting and numbering
fluid '(hcount!* pcount!* mcount!* fcount!* bcount!* b4count!*
basecount!* hzerocount!*);
% control of the polynomial arithmetic actually loaded
fluid '(currentvdpmodule!*);
fluid '(glexmat!*); % matrix for the indirect lex ordering
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% interface functions
% parameters;
% glexconvert(basis,[vars],[maxdeg=n],[newvars={x,y,..}])
symbolic procedure glexconverteval u;
begin integer n; scalar !*groebfac,!*groebrm,!*factor,!*gsugar,
v, bas,vars,maxdeg,newvars,!*exp; !*exp := t;
u := for each p in u collect reval p;
bas := car u ; u := cdr u;
while u do
<< v := car u; u := cdr u;
if eqcar(v,'list) and null vars then vars := v
else if eqcar(v,'equal) then
if(v := cdr v)and eqcar(v,'maxdeg) then maxdeg := cadr v
else if eqcar(v,'newvars) then newvars := cadr v
else <<prin2 (car v);
rerror(groebnr2,4,"Glexconvert, keyword unknown")>>
else rerror(groebnr2,5,
"Glexconvert, too many positional parameters")>>;
return glexbase1(bas,vars,maxdeg,newvars);
end;
put('glexconvert,'psopfn,'glexconverteval);
symbolic procedure glexbase1(u,v,maxdeg,nv);
begin scalar vars,w,nd,oldorder,!*gcd,!*ezgcd,!*gsugar;
integer pcount!*;
!*gcd := t;
w := for each j in groerevlist u
collect if eqexpr j then !*eqn2a j else j;
if null w then rerror(groebnr2,6,"Empty list in Groebner");
vars := groebnervars(w,v);
!*vdpinteger := !*vdpmodular := nil;
if not flagp(dmode!*,'field) then !*vdpinteger := t
else
if !*modular then !*vdpmodular := t;
if null vars then vdperr 'groebner;
oldorder := vdpinit vars;
% cancel common denominators
w := for each j in w collect reorder numr simp j;
% optimize varable sequence if desired
w := for each j in w collect f2vdp j;
for each p in w do
nd := nd or not vdpcoeffcientsfromdomain!? p;
if nd then <<!*vdpmodular:= nil;
!*vdpinteger := t;
glexdomain!* := 2>>
else glexdomain!* := 1;
if glexdomain!*=1 and not !*vdpmodular then !*ezgcd:=t;
if null maxdeg then maxdeg := 200;
if nv then nv := groerevlist nv;
if null nv then nv := vars else
for each x in nv do if not member(x,vars) then
<<rerror(groebnr2,7,list("new variable ",x,
" is not a basis variable"))>>;
u := for each v in nv collect a2vdp v;
gbtest w;
w := glexbase2(w,u,maxdeg);
w := 'list . for each j in w collect prepf j;
setkorder oldorder;
gvarslast := 'list . vars;
return w;
end;
fluid '(glexeqsys!* glexvars!* glexcount!* glexsub!*);
symbolic procedure glexbase2(oldbase,vars,maxdeg);
% in contrast to documented algorithm monbase ist a list of
% triplets (mon.cof.vect)
% such that cof * mon == vect modulo oldbase
% (cof is needed because of the integer algoritm)
begin scalar lexbase, staircase, monbase;
scalar monom, listofnexts, vect,q,glexeqsys!*,glexvars!*,
glexsub!*;
integer n;
if not groezerodim!?(oldbase,length vars) then
prin2t "####### warning: ideal is not zerodimensional ######";
% prepare matrix for the indirect lex ordering
glexmat!* := for each u in vars collect vdpevlmon u;
monbase := staircase := lexbase := nil;
monom := a2vdp 1;
listofnexts := nil;
while not(monom = nil) do
<<
if not glexmultipletest(monom,staircase) then
<< vect := glexnormalform(monom,oldbase);
q := glexlinrel(monom,vect,monbase);
if q then
<<lexbase := q . lexbase; maxdeg := nil;
staircase := monom . staircase;
>>
else
<<monbase := glexaddtomonbase(monom,vect,monbase);
n := n #+1;
if maxdeg and n#> maxdeg then
rerror(groebnr2,8,
"No univar. polynomial within degree bound");
listofnexts :=
glexinsernexts(monom,listofnexts,vars)>>;
>>;
if null listofnexts then monom := nil
else << monom := car listofnexts ;
listofnexts := cdr listofnexts >>;
>>;
return lexbase;
end;
symbolic procedure glexinsernexts(monom,l,vars);
begin scalar x;
for each v in vars do
<< x := vdpprod(monom,v);
if not vdpmember(x,l) then
<<vdpputprop(x,'factor,monom); vdpputprop(x,'monfac,v);
l := glexinsernexts1(x,l);
>>;
>>;
return l;
end;
symbolic procedure glexmultipletest(monom,staircase);
if null staircase then nil
else if vevmtest!?(vdpevlmon monom, vdpevlmon car staircase)
then t
else glexmultipletest(monom, cdr staircase);
symbolic procedure glexinsernexts1(m,l);
if null l then list m
else if glexcomp(vdpevlmon m,vdpevlmon car l) then m.l
else car l . glexinsernexts1(m,cdr l);
symbolic procedure glexcomp(ev1,ev2);
% true if ev1 is greater than ev2
% we use an indirect ordering here (mapping via newbase variables)
glexcomp0(glexcompmap(ev1,glexmat!*),
glexcompmap(ev2,glexmat!*));
symbolic procedure glexcomp0(ev1,ev2);
if null ev1 then nil
else if null ev2 then glexcomp0(ev1,'(0))
else if (car ev1 #- car ev2) = 0
then glexcomp0(cdr ev1,cdr ev2)
else if car ev1 #< car ev2 then t
else nil;
symbolic procedure glexcompmap (ev,ma);
if null ma then nil
else glexcompmap1(ev,car ma) . glexcompmap(ev,cdr ma);
symbolic procedure glexcompmap1(ev1,ev2);
% the dot product of two vectors
if null ev1 or null ev2 then 0
else (car ev1 #* car ev2) #+ glexcompmap1(cdr ev1,cdr ev2);
symbolic procedure glexaddtomonbase(monom,vect,monbase);
% primary effect: (monom . vect) . monbase
% secondary effect: builds the equation system
begin scalar x;
if null glexeqsys!* then
<<glexeqsys!* := a2vdp 0; glexcount!*:=-1>>;
x := mkid('gunivar,glexcount!*:=glexcount!*+1);
glexeqsys!* := vdpsum(glexeqsys!*,vdpprod(a2vdp x,cdr vect));
glexsub!* := (x .(monom . vect)) . glexsub!*;
glexvars!* := x . glexvars!*;
return (monom . vect) . monbase;
end;
symbolic procedure glexlinrelold(monom,vect,monbase);
if monbase then
begin scalar sys,sub,auxvars,r,v,x;
integer n;
v := cdr vect;
for each b in reverse monbase do
<<x := mkid('gunivar,n); n := n+1;
v := vdpsum(v,vdpprod(a2vdp x,cddr b));
sub := (x . b) . sub;
auxvars := x . auxvars>>;
while not vdpzero!? v do
<<sys := vdp2f vdpfmon(vdplbc v,nil) . sys; v := vdpred v>>;
x := sys;
sys := groelinsolve(sys,auxvars);
if null sys then return nil;
% construct the lex polynomial
if !*trgroeb
then prin2t " ======= constructing new basis polynomial";
r := vdp2f vdpprod(monom,car vect) ./ 1;
for each s in sub do
r:= addsq(r,multsq(vdp2f vdpprod(cadr s,caddr s) ./ 1,
cdr assoc(car s,sys)));
r := vdp2f vdpsimpcont f2vdp numr r;
return r;
end;
symbolic procedure glexlinrel(monom,vect,monbase);
if monbase then
begin scalar sys,r,v,x;
v := vdpsum(cdr vect,glexeqsys!*);
while not vdpzero!? v do
<<sys := vdp2f vdpfmon(vdplbc v,nil) . sys; v := vdpred v>>;
x := sys;
sys := groelinsolve(sys,glexvars!*);
if null sys then return nil;
% construct the lex polynomial
r := vdp2f vdpprod(monom,car vect) ./ 1;
for each s in glexsub!* do
r:= addsq(r,multsq(vdp2f vdpprod(cadr s,caddr s) ./ 1,
cdr assoc(car s,sys)));
r := vdp2f vdpsimpcont f2vdp numr r;
return r;
end;
symbolic procedure glexnormalform(m,g);
% reduce m wrt basis G;
% the reduction product is preserved in m for later usage
begin scalar cof,vect,r,f,fac1;
if !*trgroeb then prin2t " ======= reducing ";
fac1 := vdpgetprop(m,'factor);
if fac1 then vect := vdpgetprop(fac1,'vector);
if vect then
<< f := vdpprod(cdr vect,vdpgetprop(m,'monfac));
cof := car vect>>
else
<< f := m; cof := a2vdp 1>>;
r := glexnormalform1(f,g,cof);
vdpputprop(m,'vector,r);
if !*trgroeb then
<<vdpprint vdpprod(car r,m); prin2t " =====> ";
vdpprint cdr r>>;
return r;
end;
symbolic procedure glexnormalform1(f,g,cof);
begin scalar f1,c,vev,divisor,done,fold,a,b;
fold := f; f1 := vdpzero(); a:= a2vdp 1;
while not vdpzero!? f do
begin
vev:=vdpevlmon f; c:=vdplbc f;
divisor := groebsearchinlist (vev,g);
if divisor then done := t;
if divisor then
if !*vdpinteger then
<< f:=groebreduceonestepint(f,a,c,vev,divisor);
b := secondvalue!*;
cof := vdpprod(b,cof);
if not vdpzero!? f1 then
f1 := vdpprod(b,f1);
>>
else
f := groebreduceonesteprat(f,nil,c,vev,divisor)
else
<<f1 := vdpappendmon (f1,vdplbc f,vdpevlmon f);
f := vdpred f;
>>;
end;
if not done then return cof . fold;
f := groebsimpcont2(f1,cof); cof := secondvalue!*;
return cof . f;
end;
symbolic procedure groelinsolve(equations,xvars);
(begin scalar r,q,test,oldmod,oldmodulus;
if !*trgroeb then prin2t " ======= testing linear dependency ";
r := t;
if not !*modular and glexdomain!*=1 then
<<oldmod := dmode!*;
if oldmod then setdmode(get(oldmod,'dname),nil);
oldmodulus := current!-modulus;
setmod list 16381; % = 2**14-3
setdmode('modular,t);
r := groelinsolve1(for each u in equations collect
numr simp prepf u, xvars);
setdmode('modular,nil);
setmod list oldmodulus;
if oldmod then setdmode(get(oldmod,'dname),t);
>> where !*ezgcd=nil;
if null r then return nil;
r := groelinsolve1(equations,xvars);
if null r then return nil;
% divide out the common content
for each s in r do
if not(denr cdr s = 1) then test := t;
if test then return r;
q := numr cdr car r;
% for each s in cdr r do
% if q neq 1 then
% q := gcdf!*(q,numr cdr s);
% if q=1 then return r;
% r := for each s in r collect
% car s . (quotf(numr cdr s, q) ./ 1);
return r;
end) where !*ezgcd=!*ezgcd; % stack old value
symbolic procedure groelinsolve1(equations,xvars);
% Gaussian elimination in integer mode
% free of unexact divisions (see Davenport et al,CA, pp86-87
% special cases: trivial equations are ruled out early
% INPUT:
% equations: list of standard forms
% xvars: variables for the solution
% OUTPUT:
% list of pairs (var . solu) where solu is a standard quotient
%
% internal data structure: standard forms as polynomials in xvars
begin scalar oldorder,x,p,solutions,val,later,break,gc,field;
oldorder := setkorder xvars;
field := dmode!* and flagp(dmode!*,'field);
equations := for each eqa in equations collect reorder eqa;
for each eqa in equations do
if eqa and domainp eqa then break:= t;
if break then goto empty;
equations := sort(equations,function grloelinord);
again:
break := nil;
for each eqa in equations do if not break then
% car step: eliminate equations of type 23 = 0
% and 17 * u = 0
% and 17 * u + 22 = 0;
<< if null eqa then equations := delete(eqa,equations)
else if domainp eqa then break := t % inconsistent system
else if not member(mvar eqa,xvars) then break := t
else if domainp red eqa or not member(mvar red eqa,xvars) then
<<equations := delete (eqa,equations);
x := mvar eqa;
val := if lc eqa = 1 then negf red eqa ./ 1
else multsq(negf red eqa ./ 1, 1 ./lc eqa);
solutions := (x . val) . solutions;
equations := for each q in equations collect
groelinsub(q,list(x.val));
later :=
for each q in later collect groelinsub(q,list(x.val));
break := 0;
>>;
>>;
if break = 0 then goto again else if break then goto empty;
% perform an elimination loop
if null equations then goto ready;
equations := sort(equations,function grloelinord);
p := car equations; x := mvar p;
equations := for each eqa in cdr equations collect
if mvar eqa = x then
<< if field then
eqa := addf(eqa, negf multf(quotf(lc eqa,lc p),p))
else
<<gc := gcdf(lc p,lc eqa);
eqa := addf(multf(quotf(lc p,gc),eqa),
negf multf(quotf(lc eqa,gc),p))>>;
if not domainp eqa then
eqa := numr multsq(eqa ./ 1, 1 ./ lc eqa);
%%%%%%eqa := groelinscont(eqa,xvars);
eqa>>
else eqa;
later := p . later;
goto again;
ready: % do backsubstitutions
while later do
<<p := car later; later := cdr later;
p := groelinsub(p,solutions);
if domainp p or not member(mvar p,xvars) or
(not domainp red p and member(mvar red p,xvars)) then
<<break := t; later := nil>>;
x := mvar p;
val := if lc p = 1 then negf red p ./ 1
else quotsq(negf red p ./ 1 , lc p ./ 1);
solutions := (x . val) . solutions;
>>;
if break then goto empty else goto finis;
empty: solutions := nil;
finis: setkorder oldorder;
solutions := for each s in solutions collect
car s . (reorder numr cdr s ./ reorder denr cdr s);
return solutions;
end;
symbolic procedure grloelinord(u,v);
% apply ordop to the mainvars of u and v
ordop(mvar u, mvar v);
symbolic procedure groelinscont(f,vars);
% reduce content from standard form f
if domainp f then f else
begin scalar c;
c := groelinscont1(lc f,red f,vars);
if c=1 then return f;
prin2 "*************content: ";print c;
return quotf(f,c);
end;
symbolic procedure groelinscont1(q,f,vars);
% calculate the contents of standard form f
if null f or q=1 then q
else if domainp f or not member(mvar f,vars) then gcdf!*(q,f)
else groelinscont1(gcdf!*(q,lc f),red f,vars);
symbolic procedure groelinsub(s,a);
% s is a standard form linear in the top level variables
% a is an assiciation list (variable . sq) . ...
% The value is the standard form, where all substitutions
% from a are done in s (common denominator ignored).
numr groelinsub1(s,a);
symbolic procedure groelinsub1(s,a);
if domainp s then s ./ 1
else (if x then addsq(multsq(cdr x,lc s ./ 1),y)
else addsq(lt s .+ nil ./ 1,y))
where x=assoc(mvar s,a),y=groelinsub1(red s,a);
endmodule;
end;