Artifact 35219583bcaad098525d257ef9b35b7e2eb5bcc13026025d755812b11a8b6af8:
- Executable file
r37/packages/groebner/glexconv.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 17721) [annotate] [blame] [check-ins using] [more...]
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;