Artifact dac1306f62a8c9c32bf98d91831b3aa4d70e485de402d471881c98057725a7b2:
- Executable file
r36/src/invbase.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: 23174) [annotate] [blame] [check-ins using] [more...]
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;