% ----------------------------------------------------------------------
% $Id: gb.red,v 1.38 2003/05/20 08:17:38 dolzmann Exp $
% ----------------------------------------------------------------------
% Copyright (c) 1999-2003 Andreas Dolzmann and Thomas Sturm
% ----------------------------------------------------------------------
% $Log: gb.red,v $
% Revision 1.38 2003/05/20 08:17:38 dolzmann
% Moved cd_init to the beginning of cgb_interface!$. This may be neccessary for the a2s procedures.
%
% Revision 1.37 2003/05/20 07:38:26 dolzmann
% Do not load modules belongig to this package.
%
% Revision 1.36 2003/05/20 07:24:46 dolzmann
% Moved macro *_mkinterface to the right place.
%
% Revision 1.35 2003/05/19 10:21:32 dolzmann
% Fixed bugs in interface code.
% Modified interface code.
%
% Revision 1.34 2003/05/05 12:48:12 dolzmann
% Added cleanup function.
% Added reval to a2s procedures.
%
% Revision 1.33 2003/05/05 11:56:06 dolzmann
% Added interface generator.
% Removed old interfaces.
% Added interface for gb_reduce.
% Added procedure gb_gbgsys for computing a non-parametric groebner system.
% Added interface for gbgsys.
% Corrected comments.
%
% Revision 1.32 2003/04/16 09:08:15 dolzmann
% Added switch !*cgbsugar.
% Corrected and added some comments.
% Remove vdp properties after gb computations.
% Added procedure gb_reduce.
%
% Revision 1.31 1999/04/15 07:08:01 dolzmann
% Do not load rltools during compilation.
%
% Revision 1.30 1999/04/13 20:57:04 dolzmann
% Renamed switches to cgb...
% Removed !*gsugar.
% Sort the input system.
%
% Revision 1.29 1999/04/13 13:51:07 dolzmann
% gb_gb was called with three arguments instead of one argument.
%
% Revision 1.28 1999/04/11 11:30:55 dolzmann
% Added standard form interface gb_gbf. Added a procedure gb_gsys!$ for
% computing non-parametric Groebner systems.
%
% Revision 1.27 1999/04/11 09:50:40 dolzmann
% Completely rewritten the interface code for the AM.
% Moved the module vdp from dp.red into this file.
% Adappted the code to the dip_init procedure.
%
% Revision 1.26 1999/04/06 11:53:57 dolzmann
% Removed switches trgroeb, trgroebr, trgreobs, and related code.
%
% Revision 1.25 1999/04/04 16:46:19 sturm
% Added copyright and CVS fluids.
%
% Revision 1.24 1999/04/04 14:50:38 sturm
% Implemented switch tdusetorder.
%
% Revision 1.23 1999/04/04 14:09:33 sturm
% Moved dip_ilcomb and dip_ilcombr from cgb.red to dp.red.
% Created vdp_ilcomb and vdp_ilcombr for gb.red.
%
% Revision 1.22 1999/04/04 12:23:16 dolzmann
% Procedure gb_spolynomial expexts a critical pair instead of two polynomials.
%
% Revision 1.21 1999/04/03 13:38:37 sturm
% Adapted to new dip_init/dip_cleanup.
% Fixed bug in gb_strange!-reduction.
%
% Revision 1.20 1999/03/31 14:09:56 sturm
% Fixed numerous details encountered during CGB reimplementation.
%
% Revision 1.19 1999/03/30 11:29:00 dolzmann
% gb_a2s returns now a list of SF's.
% Reimplemented procedure gb_vars.
%
% Revision 1.18 1999/03/30 09:36:43 dolzmann
% Removed unused procedure gb_vdpvordopt3.
%
% Revision 1.17 1999/03/30 09:34:44 dolzmann
% Procedure gb_groebner1 binds all locally used fluids.
%
% Revision 1.16 1999/03/17 12:34:06 dolzmann
% Removed fluids !*vdpinteger, !*vdpmodular, !*grmod and their bindings.
% Use vdp_monp instead of vdp_length.
% Added the lost procedure min!# under the new name gb_min!#.
%
% Revision 1.15 1999/03/12 13:29:41 sturm
% Added a procedure gb_updbase. Currently used in the style of groebner.red
% incontrast to the Asir variant. Switch gbupdb.
%
% Revision 1.14 1999/03/12 10:58:47 dolzmann
% Introduced verbose output for the length of g.
% Use procedure ev_sdivp of package dp instead of gb_vevsdivp.
% Removed test for zero polynomials in gb_spolynomial.
% Moved procedure vdp_ilcomb1 and vdp_ilcomb1r into package dp and
% renamed the procedures accordingly.
%
% Revision 1.13 1999/03/12 08:57:59 dolzmann
% Introduced switch gbcheckg and related code.
% All statistics code is guarded by the respective switches.
% Replaced all calls of vdp_putprop and vdp_getprop by access functions
% for the properties of the dp package.
%
% Revision 1.12 1999/03/05 10:35:41 dolzmann
% Adapted to newly created dp package.
%
% Revision 1.11 1999/03/02 15:52:11 dolzmann
% Added switch gbcounthf for counting reducible H-polynomials with groebstat.
%
% Revision 1.10 1999/03/02 15:20:02 dolzmann
% Use private gb_vdpilcomb1 and gb_vdpilcomb1r. Only slightly faster.
%
% Revision 1.9 1999/02/25 16:04:56 sturm
% Added switch gbverbose: output pairs left, and frequency of heuristic
% content reduction when altered.
% Added switch gbcontred and fluid gb_mincontred!*.
% Switch groebstat works independently from trgroeb.
% Use ioto for printing statistics.
%
% Revision 1.8 1999/02/25 09:24:36 sturm
% Fixed memq to member in gb_tr2crit.
%
% Revision 1.7 1999/02/25 08:27:45 sturm
% Initialize gb_strangecount!* and spac in groebner2.
%
% Revision 1.6 1999/02/25 06:11:21 dolzmann
% Fixed a bug in gb_tr3crit. gb is as fast as grobener on p6_sc wrt.
% revgradlex and only 10 times slower as groebner wrt. lex.
%
% Revision 1.5 1999/02/24 10:48:43 sturm
% Added switch and fluid declarations and binding. In particular
% !*groebfullreduction, !*gtraverso!-sloppy, !*ezgcd, !*groebdivide.
% gb is faster than groebner on p6_sc now.
%
% Revision 1.4 1999/02/24 08:44:28 sturm
% Added gb_searchinlist. Successfully computes gb({x*y+1,x*z+1}).
%
% Revision 1.3 1999/02/24 08:36:13 sturm
% Added gb_vdpvordopt including patch code. Fixed calls to tt.
%
% Revision 1.2 1999/02/23 16:43:27 sturm
% Checked and corrected by AD. Compiles now.
%
% Revision 1.1 1999/02/23 16:41:38 sturm
% Initial check-in. Obtained from REDUCE 3.6 Groebner by selecting and
% reformatting.
%
% ----------------------------------------------------------------------
lisp <<
fluid '(gb_rcsid!* gb_copyright!*);
gb_rcsid!* := "$Id: gb.red,v 1.38 2003/05/20 08:17:38 dolzmann Exp $";
gb_copyright!* := "Copyright (c) 1999-2003 by A. Dolzmann and T. Sturm"
>>;
module gb;
load!-package 'ezgcd;
switch gbltbasis,groebopt,cgbstat,cgbfullred,cgbverbose,cgbcontred,
cgbcounthf,cgbcheckg;
fluid '(!*gltbasis !*groebopt !*cgbstat !*cgbfullred !*cgbverbose
!*cgbcontred !*cgbcounthf !*cgbcheckg !*cgbsloppy !*cgbsugar);
off1 'cgbstat;
on1 'cgbfullred;
off1 'cgbverbose;
off1 'cgbcontred;
off1 'cgbcounthf;
off1 'cgbcheckg;
!*cgbsloppy := T;
!*cgbsugar := nil; % Indicator for using sugar property of VDP's.
% This is on for gb computation and off for using
% the code from external procedures. Do not
% intermix this with !*gsugar from the groebner
% package. if !*cgbsugar is nil then the all
% sugars are treated as zero.
switch cgbupdb;
fluid '(cgb_updbcount!* cgb_updbcountp!* cgb_updbcalls!* !*cgbupdb);
off1 'cgbupdb;
fluid '(intvdpvars!*);
fluid '(!*gcd !*ezgcd !*factor !*exp dmode!* depl!* !*backtrace);
fluid '(secondvalue!* thirdvalue!*);
fluid '(dip_vars!* vdp_pcount!*);
fluid '(cgb_hcount!* cgb_hzerocount!* cgb_tr1count!* cgb_tr2count!*
cgb_tr3count!* cgb_b4count!* cgb_strangecount!* cgb_paircount!*
cgb_hfaccount!* cgb_gcount!* cgb_gstat!* cgb_gbcount!*);
global '(gvarslast gltb);
gltb := {'list};
fluid '(cgb_contcount!* cgb_mincontred!*);
cgb_mincontred!* := 20; % originally 10
% Generate the following interfaces for AM:
% - gb: List of Polyomials -> List of Polynomials
% - gbgsys: List of Polyomials -> Groebner System
% - reduce: Polynomial X List of Polynomials -> Polynomial
% - ltb: List of Polynomials -> List of Terms (as Polynomials)
% and the following interfaces for the SM:
% - gbf,gsysf,reducef,ltbf. All polynomials are represented as SF's,
% except the return valute of reducef, which is an SQ. All Procedures
% expect three dditional paramters: List of variable, term order, and
% additional arguments to term order.
macro procedure gb_mkinterface(argl);
begin
scalar a2sl1,a2sl2,defl,xvfn,s2a,s2s,s,modes,
args,bname,len,sm,prgn,ami,smi,psval,postfix;
bname := eval nth(argl,2);
a2sl1 := eval nth(argl,3);
a2sl2 := eval nth(argl,4);
defl := eval nth(argl,5);
xvfn := eval nth(argl,6);
s2a := eval nth(argl,7);
s2s := eval nth(argl,8);
s := eval nth(argl,9);
postfix := eval nth(argl,10);
modes := eval nth(argl,11);
len := length a2sl1;
args := for i := 1:len+3 collect mkid('a,i);
sm := intern compress append('(!g !b !_),explode bname);
% Define the symbolic mode interface
if (null modes or modes eq 'sm) then <<
smi := intern compress nconc(explode sm,explode postfix);
prgn := {'put,mkquote smi,''number!-of!-args,len+3} . prgn;
prgn := {'de,smi,args,{'gb_interface!$,mkquote sm, mkquote a2sl1,
mkquote a2sl2,mkquote defl,mkquote xvfn,mkquote s2a,mkquote s2s,
mkquote s,T,'list . args}} . prgn
>>;
if (null modes or modes eq 'am) then <<
% Define the algebraic mode interface
ami := bname;
% ami := intern compress append('(!g !b),explode bname);
psval := intern compress nconc(explode ami,'(!! !$));
prgn := {'put,mkquote ami,''psopfn,mkquote psval} . prgn;
prgn := {'put,mkquote psval,''number!-of!-args,1} . prgn;
prgn := {'put,mkquote psval,''cleanupfn,''gb_cleanup} . prgn;
prgn := {'de,psval,'(argl),{'gb_interface!$,mkquote sm, mkquote a2sl1,
mkquote a2sl2,mkquote defl,mkquote xvfn,mkquote s2a, mkquote s2s,
mkquote s,nil,'argl}} . prgn
>>;
return 'progn . prgn
end;
gb_mkinterface('gb,'(gb_a2s!-psys),'(gb_a2s2!-psys),
nil,'gb_xvars!-psys,'gb_s2a!-gbx,'gb_s2s!-gb,T,'f,nil);
gb_mkinterface('gbgsys,'(gb_a2s!-psys),'(gb_a2s2!-psys),
nil,'gb_xvars!-psys,'gb_s2a!-gsys,'gb_s2s!-gsys,T,'f,nil);
gb_mkinterface('reduce,'(gb_a2s!-pol gb_a2s!-psys),'(gb_a2s2!-pol gb_a2s2!-psys),
nil,'gb_xvars!-ppsys,'gb_s2a!-pol,'gb_s2s!-pol,T,'f,nil);
gb_mkinterface('ltb,'(gb_a2s!-psys),'(gb_a2s2!-psys),
nil,'gb_xvars!-psys,'gb_s2a!-gb,'gb_s2s!-gb,T,'f,nil);
procedure gb_a2s!-psys(l);
% Groebner bases algebraic mode to symbolic mode polynomial system.
% [l] is an AMPSYS. Returns an FPSYS.
begin scalar w,resl;
for each j in getrlist reval l do <<
w := numr simp j;
if w and not(w member resl) then
resl := w . resl
>>;
return sort(resl,'ordp)
end;
procedure gb_a2s2!-psys(fl);
for each x in fl collect vdp_f2vdp x;
procedure gb_s2a!-gb(u);
% Groebner bases symbolic mode to algebraic mode GB. [u] is a list
% of CGP's. Returns an AMPSYS.
'list . for each x in u collect vdp_2a x;
procedure gb_s2a!-gbx(l);
% symbolic to algebraic mode groebner base extended version. [l] is
% a list of VDP's. Returns an AM object. This procedure sets as a
% side effet the global variable gltb provided !*gltbasis is on.
<<
if !*gltbasis then
gltb := gb_gb2gltb l;
gb_s2a!-gb l
>>;
procedure gb_s2s!-gb(l);
gb_gb!-sfl l;
procedure gb_s2a!-gsys(l);
'list . for each x in l collect
{'list,rl_mk!*fof rl_smkn('and,car x),gb_s2a!-gb cadr x};
procedure gb_s2s!-gsys(l);
for each x in l collect
{rl_smkn('and,car x),gb_s2s!-gb cadr x};
procedure gb_a2s!-pol(p);
numr simp reval p;
procedure gb_a2s2!-pol(p);
vdp_f2vdp p;
procedure gb_s2a!-pol(p);
vdp_2a p;
procedure gb_s2s!-pol(p);
vdp_2sq p;
procedure gb_xvars!-psys(l,vl);
gb_vars(l,vl);
procedure gb_xvars!-ppsys(p,l,vl);
gb_vars(p . l,vl);
procedure gb_cleanup(u,v);
u;
procedure gb_interface!$(fname,a2sl1,a2sl2,defl,xvfn,s2a,s2s,s,smp,argl);
% fname is a function, the name of the procedure to be called;
% [a2sl1] and [as2sl2] are a list of functions, called to be
% transform algebraic arguments to symbolic arguments; [defl] is a
% list of algebraic defualt arguments; xvfn is a procedure for
% extracting the variables from all arguments; [s2a] is procedure
% for transforming the symbolic return value to an algebraic mode
% return value; [argl] is the list of arguments; [s] is a flag;
% [smp] is a flag. Return an S-expr. If [s] is on then second stage
% of argument processing is done with the results of the first one.
begin scalar w,vl,nargl,oenv,m,c,x;
if not smp then <<
nargl := gb_am!-pargl(fname,a2sl1,argl,defl);
vl := apply(xvfn,append(nargl,{td_vars()}));
oenv := vdp_init(car vl,td_sortmode(),td_sortextension());
gvarslast := 'list . car vl;
>> else <<
w := gb_sm!-pargl(argl);
nargl := car w;
m := cadr w;
c := caddr w;
x := cadddr w;
vl := apply(xvfn,append(nargl,{m}));
oenv := vdp_init(car vl,c,x);
>>;
w := errorset({'gb_interface1!$,
mkquote fname,mkquote a2sl2,mkquote s2a,mkquote s2s,mkquote s,
mkquote smp,mkquote argl, mkquote nargl,mkquote car vl,
mkquote cdr vl},T,!*backtrace);
vdp_cleanup oenv;
if errorp w then
rederr {"Error during ",fname};
return car w
end;
procedure gb_sm!-pargl(argl);
begin scalar nargl,m,c,x;
nargl := reverse argl;
x := car nargl;
nargl := cdr nargl;
c := car nargl;
nargl := cdr nargl;
m := car nargl;
nargl := cdr nargl;
nargl := reversip nargl;
return {nargl,m,c,x}
end;
procedure gb_am!-pargl(fname,a2sl1,argl,defl);
% process argument list for algebraic mode.
begin integer l1,l2,l3; scalar w,nargl,scargl,scdefl;
l1 := length argl;
l2 := length a2sl1;
l3 := l2 - length defl;
if l1 < l3 or l1 > l2 then
rederr {fname,"called with",l1,"arguments instead of",l3,"-",l2};
scargl := argl;
scdefl := defl;
nargl := for each x in a2sl1 collect <<
if scargl then <<
w := car scargl;
scargl := cdr scargl
>> else <<
w := car scdefl;
scdefl := cdr scdefl
>>;
apply(x,{w})
>>;
return nargl
end;
procedure gb_interface1!$(fname,a2sl2,s2a,s2s,s,smp,argl,nargl,m,p);
begin scalar w,pl;
pl := if s then nargl else argl;
argl := for each x in a2sl2 collect <<
w := car pl;
pl := cdr pl;
apply(x,{w})
>>;
% w := apply(fname,nconc(argl,{m,p}));
w := apply(fname,argl);
w := if smp then
apply(s2s,{w})
else
apply(s2a,{w});
return w
end;
smacro procedure gb_tt(s1,s2);
% lcm of leading terms of s1 and s2
ev_lcm(vdp_evlmon s1,vdp_evlmon s2);
procedure gb_gb!-sfl(u);
% Groebner bases GB to SF list. [u] is a list of CGP's. Returns a
% list of SF's.
for each p in u collect vdp_2f p;
procedure gb_domainchk();
% Groebner bases domain check. No argument. Return value not
% defined. Raises an error if the current domain is not valid for
% GB computations.
if not memq(dmode!*,'(nil)) then
rederr bldmsg("gb does not support domain: %w",get(dmode!*,'dname));
procedure gb_vars(l,vl); %DROPPED: depend,rules,zero divisors.
% Groebner bases variables. [l] is a list of SF's; [vl] is the list
% of main variables. Returns a pair $(m . p)$ where $m$ and $p$ are
% list of variables. $m$ is the list of used main variables and $p$
% is the list of used parameters.
begin scalar w,m,p;
for each f in l do
w := union(w,kernels f);
if vl then <<
m := gb_intersection(vl,w);
p := setdiff(w,vl)
>> else
m := w;
return gb_varsopt(l,m) . p
end;
procedure gb_intersection(a,b);
% Groebner bases intersection. [a] and [b] are lists. Returns a
% list. The returned list contains all elements occuring in [a] and
% in [b]. The order of the elements is the same as in [a].
for each x in a join
if x member b then
{x};
procedure gb_varsopt(l,vl);
% Groebner bases variables optimize. [l] is a list of SF's; [vl] is
% the list of main variables. Returns a possibly reorderd list of
% main variables.
if !*groebopt and td_sortmode() memq '(lex gradlex revgradlex) then
gb_vdpvordopt(l,vl)
else
vl;
procedure gb_gb2gltb(base);
'list . for each j in base collect
vdp_2a vdp_fmon(bc_a2bc 1,vdp_evlmon j);
procedure gb_ltb(l);
for each p in l collect
vdp_fmon(bc_a2bc 1,vdp_evlmon p);
procedure gb_gbgsys(p);
{{nil,gb_gb p,nil}};
procedure gb_gb(p);
begin scalar spac,p1,savetime,!*factor,!*exp,intvdpvars!*,!*gcd,!*ezgcd,
dip_vars!*,secondvalue!*,thirdvalue!*,cgb_gstat!*,!*cgbsugar;
integer vdp_pcount!*,cgb_contcount!*,cgb_hcount!*,cgb_hzerocount!*,
cgb_tr1count!*,cgb_tr2count!*,cgb_tr3count!*,cgb_b4count!*,
cgb_strangecount!*,cgb_paircount!*,cgb_hfaccount!*,cgb_gcount!*,
cgb_gbcount!*,cgb_updbcount!*,cgb_updbcountp!*,cgb_updbcalls!*;
!*exp := !*gcd := !*ezgcd := T;
!*cgbsugar := T;
if !*cgbstat then
savetime := time();
if !*cgbcheckg then
cgb_gstat!* := nil;
cgb_contcount!* := cgb_mincontred!*;
if !*cgbstat then
spac := gctime();
p1 := if !*cgbupdb then
gb_traverso!-sturm!-experimental p
else
gb_traverso p;
if !*cgbstat then <<
ioto_tprin2t "Statistics for GB computation:";
ioto_prin2t {"Time: ",time() - savetime," ms plus GC time: ",
gctime() - spac," ms"};
ioto_prin2t {"H-polynomials total: ",cgb_hcount!*};
ioto_prin2t {"H-polynomials zero: ",cgb_hzerocount!*};
if !*cgbcounthf then
ioto_prin2t {"H-polynomials reducible: ",cgb_hfaccount!*};
if !*cgbcheckg then
ioto_prin2t {"H-polynomials gaussible: ",cgb_gcount!*,
" ",cgb_gstat!*};
ioto_prin2t {"Crit Tr1 hits: ",cgb_tr1count!*};
ioto_prin2t {"Crit B4 hits: ",cgb_b4count!*," (Buchberger 1)"};
ioto_prin2t {"Crit Tr2 hits: ",cgb_tr2count!*};
ioto_prin2t {"Crit Tr3 hits: ",cgb_tr3count!*};
if !*cgbupdb then
ioto_prin2t {"updbase: calls ",cgb_updbcalls!*,", del ",
cgb_updbcountp!*,"/",cgb_updbcount!*};
ioto_prin2t {"Strange reductions: ",cgb_strangecount!*}
>>;
return p1
end;
procedure gb_traverso!-sturm!-experimental(g0);
begin scalar gall,g,d,s,h,p;
g0 := for each fj in g0 join
if not vdp_zero!? fj then
{vdp_setsugar(vdp_enumerate vdp_simpcont fj,vdp_tdeg fj)};
for each h in g0 do << % create initial critical pairs
p := {nil,h,h};
h := gb_enumerate h;
d := gb_traverso!-pairlist(h,g,d);
g := nconc(g,{h});
gall := nconc(gall,{h});
>>;
while d do << % critical pairs left
if !*cgbverbose then <<
ioto_prin2 {"[",cgb_paircount!*,"] "};
cgb_paircount!* := cgb_paircount!* #- 1
>>;
p := car d;
d := cdr d;
s := gb_spolynomial(p);
h := gb_simpcontnormalform gb_normalform(s,gall);
cgb_hcount!* := cgb_hcount!* #+ 1;
if vdp_zero!? h then <<
cgb_hzerocount!* := cgb_hzerocount!* #+ 1;
>> else if ev_zero!? vdp_evlmon h then << % base 1 found
h := gb_enumerate h;
d := nil;
g := {h}
>> else <<
if !*cgbcounthf then
if cddr fctrf vdp_2f h then
cgb_hfaccount!* := cgb_hfaccount!* #+ 1;
h := gb_enumerate h;
d := gb_traverso!-pairlist(h,g,d);
gall := gb_updbase(gall,h);
g := nconc(g,{h})
>>
>>;
return gb_traverso!-final g
end;
procedure gb_traverso(g0);
begin scalar g,d,s,h,p,gstat;
g0 := for each fj in g0 join
if not vdp_zero!? fj then
{vdp_setsugar(vdp_enumerate vdp_simpcont fj,vdp_tdeg fj)};
for each h in g0 do << % create initial critical pairs
p := {nil,h,h};
h := gb_enumerate h;
d := gb_traverso!-pairlist(h,g,d);
g := nconc(g,{h})
>>;
if !*cgbverbose then
cgb_gbcount!* := length g;
while d do << % critical pairs left
if !*cgbverbose then <<
ioto_prin2 {"[",cgb_paircount!*,"] "};
cgb_paircount!* := cgb_paircount!* #- 1
>>;
p := car d;
d := cdr d;
s := gb_spolynomial(p);
h := gb_simpcontnormalform gb_normalform(s,g);
if !*cgbstat then
cgb_hcount!* := cgb_hcount!* #+ 1;
if vdp_zero!? h then
(if !*cgbstat then
cgb_hzerocount!* := cgb_hzerocount!* #+ 1)
else if vdp_unit!? h then <<
h := gb_enumerate h;
d := nil;
g := {h}
>> else <<
if !*cgbcounthf then
if cddr fctrf vdp_2f h then
cgb_hfaccount!* := cgb_hfaccount!* #+ 1;
if !*cgbcheckg then <<
gstat := gb_chkgauss vdp_poly h;
if 1 member gstat then <<
cgb_gcount!* := cgb_gcount!* #+ 1;
cgb_gstat!* := gb_chkgauss!-stat2vl gstat . cgb_gstat!*
>>
>>;
h := gb_enumerate h;
d := gb_traverso!-pairlist(h,g,d);
g := nconc(g,{h});
if !*cgbverbose then
cgb_gbcount!* := cgb_gbcount!* #+ 1
>>
>>;
return gb_traverso!-final g
end;
procedure gb_updbase(g,h);
begin scalar hev,oc;
if !*cgbstat then
oc := cgb_updbcountp!*;
hev := vdp_evlmon h;
g := for each p in g join
if not ev_divides!?(hev,vdp_evlmon p) then
{p}
else <<
if !*cgbverbose then
ioto_prin2 "#";
if !*cgbstat then
cgb_updbcountp!* := cgb_updbcountp!* #+ 1;
nil
>>;
if !*cgbstat then <<
if not (oc #= cgb_updbcountp!*) then
cgb_updbcount!* := cgb_updbcount!* #+ 1;
cgb_updbcalls!* := cgb_updbcalls!* + 1
>>;
return nconc(g,{h})
end;
procedure gb_chkgauss(p);
% [p] is a dipoly.
begin scalar stat;
stat := for each x in dip_vars!* collect 0; %TODO: Reference to global var
while p do <<
stat := gb_chkgauss1(dip_evlmon p,stat);
p := dip_mred p
>>;
return stat
end;
procedure gb_chkgauss1(ev,stat);
begin scalar nstat,e,s,td;
td := ev_tdeg ev;
if td = 0 then
return stat;
td := if td #= 1 then 1 else -1;
while ev do <<
e := car ev;
ev := cdr ev;
s := car stat;
stat := cdr stat;
if e #=0 then
nstat := s . nstat
else if e #> 1 or s #= -1 then
nstat := (-1) . nstat
else if e #= 1 and s #= 0 then
nstat := td . nstat
else if s #=0 and e #=0 then
nstat := 0 . nstat
else
rederr "ich sehe es anders"
>>;
return reversip nstat
end;
procedure gb_chkgauss!-stat2vl(gstat);
begin scalar scdv,r;
scdv := dip_vars!*;
while gstat do <<
if eqcar(gstat,1) then
r := car scdv . r;
gstat := cdr gstat;
scdv := cdr scdv
>>;
return reversip r
end;
procedure gb_enumerate(f);
% f is a temporary result. Prepare it for medium range storage, and
% assign a number.
if vdp_zero!? f then
f
else <<
vdp_condense f;
if vdp_number f #= 0 then
f := vdp_setnumber(f,vdp_pcount!* := vdp_pcount!* #+ 1);
f
>>;
procedure gb_traverso!-pairlist(gk,g,d);
% gk: new polynomial, g: current basis, d: old pair list.
begin scalar ev,r,n;
d := gb_traverso!-pairs!-discard1(gk,d);
% build new pair list:
ev := vdp_evlmon gk;
for each p in g do
if not gb_buchcrit4t(ev,vdp_evlmon p) then <<
if !*cgbstat then
cgb_b4count!* := cgb_b4count!* #+ 1;
r := ev_lcm(ev,vdp_evlmon p) . r
>> else
n := gb_makepair(p,gk) . n;
n := gb_tr2crit(n,r);
n := gb_cplistsort(n,!*cgbsloppy);
n := gb_tr3crit n;
if !*cgbverbose and n then <<
cgb_paircount!* := cgb_paircount!* #+ length n;
ioto_cterpri();
ioto_prin2 {"(",cgb_gbcount!*,") "}
>>;
return gb_cplistmerge(d,reversip n)
end;
procedure gb_tr2crit(n,r);
% delete equivalents to coprime lcm
for each p in n join
if ev_member(car p,r) then <<
if !*cgbstat then
cgb_tr2count!* := cgb_tr2count!* #+ 1;
nil
>> else
{p};
procedure gb_tr3crit(n);
begin scalar newn,scannewn,q;
for each p in n do <<
scannewn := newn;
q := nil;
while scannewn do
if ev_divides!?(caar scannewn,car p) then <<
q := t;
scannewn := nil;
if !*cgbstat then
cgb_tr3count!* := cgb_tr3count!* #+ 1
>> else
scannewn := cdr scannewn;
if not q then
newn := gb_cplistsortin(p,newn,nil)
>>;
return newn
end;
procedure gb_traverso!-pairs!-discard1(gk,d);
% crit B. Delete triange relations.
for each pij in d join
if gb_traverso!-trianglep(cadr pij,caddr pij,gk,car pij) then <<
if !*cgbstat then
cgb_tr1count!* := cgb_tr1count!* #+ 1;
if !*cgbverbose then
cgb_paircount!* := cgb_paircount!* #- 1;
nil
>> else
{pij};
procedure gb_traverso!-trianglep(gi,gj,gk,tij);
ev_sdivp(gb_tt(gi,gk),tij) and ev_sdivp(gb_tt(gj,gk),tij);
procedure gb_traverso!-final(g);
% Final reduction and sorting.
for each rg on vdp_lsort g join
if not gb_searchinlist(vdp_evlmon car rg,cdr rg) then
{vdp_remplist gb_simpcontnormalform gb_normalform(car rg,cdr rg)};
procedure gb_buchcrit4t(e1,e2);
% nonconstructive test of lcm(e1,e2) = e1 + e2 equivalent: no
% matches of nonzero elements.
not ev_disjointp(e1,e2);
procedure gb_cplistsort(g,sloppy);
begin scalar gg;
for each p in g do
gg := gb_cplistsortin(p,gg,sloppy);
return gg
end;
procedure gb_cplistsortin(p,pl,sloppy);
% Distributive polynomial critical pair list sort. pl is a special
% list for Groebner calculation, p is a pair. Returns the updated
% list pl (p sorted into).
if null pl then
{p}
else <<
gb_cplistsortin1(p,pl,sloppy);
pl
>>;
procedure gb_cplistsortin1(p,pl,sloppy);
% Destructive insert of p into nonnull pl.
if not gb_cpcompless!?(car pl,p,sloppy) then <<
rplacd(pl,car pl . cdr pl);
rplaca(pl,p)
>> else if null cdr pl then
rplacd(pl,{p})
else
gb_cplistsortin1(p,cdr pl,sloppy);
procedure gb_cpcompless!?(p1,p2,sloppy);
% Compare 2 pairs wrt. their sugar (=cadddr) or their lcm (=car).
if sloppy then
ev_compless!?(car p1,car p2)
else
gb_cpcompless!?s(p1,p2,cadddr p1 #- cadddr p2,ev_comp(car p1,car p2));
procedure gb_cpcompless!?s(p1,p2,d,q);
if d neq 0 then
d #< 0
else if q neq 0 then
q < 0
else
vdp_number(caddr p1) #< vdp_number(caddr p2);
procedure gb_cplistmerge(pl1,pl2);
% Distributive polynomial critical pair list merge. pl1 and pl2 are
% critical pair lists used in the Groebner calculation.
% groebcplistmerge(pl1,pl2) returns the merged list.
begin scalar cpl1,cpl2;
if null pl1 then
return pl2;
if null pl2 then
return pl1;
cpl1 := car pl1;
cpl2 := car pl2;
return if gb_cpcompless!?(cpl1,cpl2,nil) then
cpl1 . gb_cplistmerge(cdr pl1,pl2)
else
cpl2 . gb_cplistmerge(pl1,cdr pl2)
end;
procedure gb_makepair(f,h);
% Construct a pair from polynomials f and h.
begin scalar ttt,sf,sh;
ttt := gb_tt(f,h);
sf := vdp_sugar(f) #+ ev_tdeg ev_dif(ttt,vdp_evlmon f);
sh := vdp_sugar(h) #+ ev_tdeg ev_dif(ttt,vdp_evlmon h);
return {ttt,f,h,ev_max!#(sf,sh)}
end;
procedure gb_spolynomial(pr);
begin scalar p1,p2,s;
p1 := cadr pr;
p2 := caddr pr;
s := gb_spolynomial1(p1,p2); % TODO: Switch for strange reduction
if vdp_zero!? s or vdp_unit!? s then
return s;
% return vdp_setsugar(gb_strange!-reduction(s,p1,p2),cadddr pr)
return gb_strange!-reduction(s,p1,p2) % TODO: normal suger for
% special cases.
end;
procedure gb_spolynomial1(p1,p2);
begin scalar ep1,ep2,ep,rp1,rp2,db1,db2,x;
ep1 := vdp_evlmon p1;
ep2 := vdp_evlmon p2;
ep := ev_lcm(ep1,ep2);
rp1 := vdp_mred p1;
rp2 := vdp_mred p2;
if vdp_zero!? rp1 and vdp_zero!? rp2 then
return rp1;
if vdp_zero!? rp1 then
return vdp_prod(rp2,vdp_fmon(bc_a2bc 1,ev_dif(ep,ep2)));
if vdp_zero!? rp2 then
return vdp_prod(rp1,vdp_fmon(bc_a2bc 1,ev_dif(ep,ep1)));
db1 := vdp_lbc p1;
db2 := vdp_lbc p2;
x := bc_gcd(db1,db2);
if not bc_one!? x then <<
db1 := bc_quot(db1,x);
db2 := bc_quot(db2,x)
>>;
return vdp_ilcomb(rp2,db1,ev_dif(ep,ep2),rp1,bc_neg db2,ev_dif(ep,ep1))
end;
procedure gb_strange!-reduction(s,p1,p2);
% Subtracts multiples of p2 from the s-polynomial such that head
% terms are eliminated early.
begin scalar tp1,tp2,ts,c,saves;
saves := s;
tp1 := vdp_evlmon p1;
tp2 := vdp_evlmon p2;
c := T; while c and not vdp_zero!? s do <<
ts := vdp_evlmon s;
if gb_buch!-ev_divides!?(tp2,ts) then
s := gb_reduceonestepint(s,vdp_zero(),vdp_lbc s,ts,p2)
else if gb_buch!-ev_divides!?(tp1,ts) then
s := gb_reduceonestepint(s,vdp_zero(),vdp_lbc s,ts,p1)
else
c := nil
>>;
if !*cgbstat and not (s eq saves) then
cgb_strangecount!* := cgb_strangecount!* #+ 1;
return s
end;
procedure gb_buch!-ev_divides!?(vev1,vev2);
% Test if vev1 divides vev2 for exponent vectors vev1 and vev2.
ev_mtest!?(vev2,vev1);
procedure gb_normalform(f,g);
% General procedure for the reduction of one polynomial modulo a
% set. f is a polynomial, g is a set of polynomials. Procedure
% behaves like type='list. f is reduced modulo g. f is possibly
% multiplied by an factor.
begin scalar f1,c,vev,divisor,tai,fold; integer n;
fold := f;
f1 := vdp_setsugar(vdp_zero(),vdp_sugar f);
while not vdp_zero!? f do <<
vev := vdp_evlmon f;
c := vdp_lbc f;
divisor := gb_searchinlist(vev,g);
if divisor then <<
tai := T;
if vdp_monp divisor then
f := vdp_cancelmev(f,vdp_evlmon divisor)
else <<
f := gb_reduceonestepint(f,f1,c,vev,divisor);
f1 := secondvalue!*;
if !*cgbcontred then <<
f := gb_adtssimpcont(f,f1,n);
f1 := secondvalue!*;
n := thirdvalue!*
>>
>>
>> else if !*cgbfullred then <<
f := gb_shift(f,f1);
f1 := secondvalue!*
>> else <<
f1 := vdp_sum(f1,f);
f := vdp_zero()
>>
>>;
return if tai then f1 else fold
end;
procedure gb_searchinlist(vev,g);
% search for a polynomial in the list G, such that the lcm divides
% vev; G is expected to be sorted in descending sequence.
if null g then
nil
else if gb_buch!-ev_divides!?(vdp_evlmon car g, vev) then
car g
else
gb_searchinlist(vev,cdr g);
procedure gb_adtssimpcont(f,f1,n);
begin scalar f0;
if vdp_zero!? f then <<
secondvalue!* := f1;
thirdvalue!* := 0;
return f
>>;
n := n + 1;
if n #> cgb_contcount!* then <<
f0 := f;
f := gb_simpcont2(f,f1);
f1 := secondvalue!*;
gb_contentcontrol(f neq f0);
n := 0
>>;
secondvalue!* := f1;
thirdvalue!* := n;
return f
end;
procedure gb_simpcont2(f,f1);
% Simplify two polynomials with the gcd of their contents. [f] is
% not zero.
begin scalar c,s1,s2;
c := vdp_content f;
if bc_one!? bc_abs c then <<
secondvalue!* := f1;
return f
>>;
s1 := vdp_sugar f;
s2 := vdp_sugar f1;
if not vdp_zero!? f1 then <<
c := vdp_content1(f1,c);
if bc_one!? bc_abs c then <<
secondvalue!* := f1;
return f
>>;
f1 := vdp_bcquot(f1,c)
>>;
f := vdp_bcquot(f,c);
vdp_setsugar(f,s1);
vdp_setsugar(f1,s2);
secondvalue!* := f1;
return f
end;
procedure gb_contentcontrol(u);
% u indicates, that a substantial content reduction was done;
% update content reduction limit from u.
<<
cgb_contcount!* := if u then
ev_max!#(0,cgb_contcount!* #- 1)
else
gb_min!#(cgb_mincontred!*,cgb_contcount!* #+ 1);
if !*cgbverbose then
ioto_prin2 {"<",cgb_contcount!*,"> "}
>>;
procedure gb_min!#(a,b);
if a #< b then a else b;
procedure gb_shift(f,f1);
begin scalar s,s1;
s := vdp_sugar f;
s1 := vdp_sugar f1;
f1 := vdp_nconcmon(f1,vdp_lbc f,vdp_evlmon f);
f := vdp_mred f;
vdp_setsugar(f,s);
vdp_setsugar(f1,ev_max!#(s1,s));
secondvalue!* := f1;
return f
end;
procedure gb_reduceonestepint(f,f1,c,vev,g1);
% Reduction step for integer case: calculate f = a*f - b*g a, b
% such that leading term vanishes (vev of lvbc g divides vev of
% lvbc f) and calculate f1 = a*f1. Return value=f, secondvalue=f1.
begin scalar vevcof,a,b,cg,x,rg1;
rg1 := vdp_mred g1;
if vdp_zero!? rg1 then << % g1 is monomial
f := vdp_mred f;
secondvalue!* := f1;
return f
>>;
vevcof := ev_dif(vev,vdp_evlmon g1); % nix lcm
cg := vdp_lbc g1;
x := bc_gcd(c,cg);
a := bc_quot(cg,x);
b := bc_quot(c,x);
% multiply relevant parts of f and f1 by a (vbc)
if not vdp_zero!? f1 then
f1 := vdp_bcprod(f1,a);
f := vdp_ilcombr(vdp_mred f,a,rg1,bc_neg b,vevcof);
secondvalue!*:= f1;
return f
end;
procedure gb_simpcontnormalform(h);
% simpCont version preserving the property SUGAR.
if vdp_zero!? h then
h
else
vdp_setsugar(vdp_simpcont h,vdp_sugar h);
procedure gb_vdpvordopt(w,vars);
% w : list of polynomials (standard forms)
% vars: list of variables
% return vars
begin scalar c;
vars := sort(vars,'ordop);
c := for each x in vars collect x . 0 . 0;
for each poly in w do
gb_vdpvordopt1(poly,vars,c);
c := sort(c,function gb_vdpvordopt2);
intvdpvars!* := for each v in c collect car v;
vars := gb_vdpvordopt31 intvdpvars!*;
return vars
end;
procedure gb_vdpvordopt31(u);
begin scalar v,y;
if null u then
return nil;
v := for each x in u join <<
y := assoc(x,depl!*);
if null y or null xnp(cdr y,u) then
{x}
>>;
return nconc(gb_vdpvordopt31 setdiff(u,v),v)
end;
procedure gb_vdpvordopt1(p,vl,c);
if null p then
0
else if domainp p or null vl then
1
else if mvar p neq car vl then
gb_vdpvordopt1(p,cdr vl,c)
else begin scalar var,pow,slot; integer n;
n := gb_vdpvordopt1 (lc p,cdr vl,c);
var := mvar p;
pow := ldeg p;
slot := assoc(var,c);
if pow #> cadr slot then <<
rplaca(cdr slot,pow);
rplacd(cdr slot,n)
>> else
rplacd(cdr slot,n #+ cddr slot);
return n #+ gb_vdpvordopt1 (red p,vl,c)
end;
procedure gb_vdpvordopt2(sl1,sl2);
% compare two slots from the power table
<<
sl1 := cdr sl1;
sl2 := cdr sl2;
car sl1 #< car sl2 or car sl1 = car sl2 and cdr sl1 #< cdr sl2
>>;
endmodule; % gb
module vdp;
%DS
% VDP ::= ('vdp,EV,BC,DP,PLIST)
%DS TERM
% A VDP, such that the polynomial contains one monomial with base
% coefficient 1.
%DS MONOMIAL
% A VDP, such that the polynomial contains one monomial.
procedure vdp_lbc(u);
caddr u;
procedure vdp_evlmon(u);
cadr u;
procedure vdp_poly(u);
car cdddr u;
procedure vdp_zero!?(u);
null vdp_poly u;
procedure vdp_plist(u);
cadr cdddr u;
procedure vdp_remplist(u);
% virtual distributive polynomial remove plist. [u] is a VDP.
% Returns a VDP. Sets the plist of [u] in-plce to [nil].
<<
nth(u,5) := nil;
u
>>;
procedure vdp_number(f);
vdp_getprop(f,'number) or 0;
procedure vdp_sugar(f);
if (vdp_zero!? f or not(!*cgbsugar)) then 0 else vdp_getprop(f,'sugar) or 0;
procedure vdp_unit!?(p);
not vdp_zero!? p and ev_zero!? vdp_evlmon p;
procedure vdp_make(vbc,vev,form);
{'vdp,vev,vbc,form,nil};
procedure vdp_monp(u);
dip_monp vdp_poly u;
procedure vdp_tdeg(u);
dip_tdeg vdp_poly u;
procedure vdp_fdip(u);
if null u then
vdp_zero()
else
vdp_make(dip_lbc u,dip_evlmon u,u);
procedure vdp_appendmon(vdp,coef,vev);
% Add a monomial to the end of a vdp (vev remains unchanged).
if vdp_zero!? vdp then
vdp_fmon(coef,vev)
else if bc_zero!? coef then
vdp
else
vdp_make(vdp_lbc vdp,vdp_evlmon vdp,dip_appendmon(vdp_poly vdp,coef,vev));
procedure vdp_nconcmon(vdp,coef,vev);
% Add a monomial to the end of a vdp (vev remains unchanged).
if vdp_zero!? vdp then
vdp_fmon(coef,vev)
else if bc_zero!? coef then
vdp
else
vdp_make(vdp_lbc vdp,vdp_evlmon vdp,dip_nconcmon(vdp_poly vdp,coef,vev));
procedure vdp_bcquot(p,c);
begin scalar r;
r := vdp_fdip dip_bcquot(vdp_poly p,c);
vdp_setsugar(r,vdp_sugar p);
return r
end;
procedure vdp_content(p);
dip_contenti vdp_poly p;
procedure vdp_content1(d,c);
dip_contenti1(vdp_poly d,c);
procedure vdp_length(f);
dip_length vdp_poly f;
procedure vdp_bcprod(p,b);
begin scalar r;
r := vdp_fdip dip_bcprod(vdp_poly p,b);
vdp_setsugar(r,vdp_sugar p);
return r
end;
procedure vdp_cancelmev(p,vev);
begin scalar r;
r := vdp_fdip dip_cancelmev(vdp_poly p,vev);
vdp_setsugar(r,vdp_sugar p);
return r
end;
procedure vdp_sum(d1,d2);
begin scalar r;
r := vdp_fdip dip_sum(vdp_poly d1,vdp_poly d2);
vdp_setsugar(r,ev_max!#(vdp_sugar d1,vdp_sugar d2));
return r
end;
procedure vdp_prod(d1,d2);
begin scalar r;
r := vdp_fdip dip_prod(vdp_poly d1,vdp_poly d2);
vdp_setsugar(r,vdp_sugar d1 #+ vdp_sugar d2);
return r
end;
procedure vdp_zero();
vdp_make('invalid,'invalid,nil);
procedure vdp_mred(u);
begin scalar r;
r := dip_mred vdp_poly u;
if null r then
return vdp_zero();
r := vdp_make(dip_lbc r,dip_evlmon r,r);
vdp_setsugar(r,vdp_sugar u);
return r
end;
procedure vdp_condense(f);
dip_condense vdp_poly f;
procedure vdp_setsugar(p,s);
% virtual distributive polynomial set sugar. [p] is a VDP, s is a
% machine integer. Returns a VDP. The sugar
% property of [p] is set to [s].
if not !*cgbsugar then p else vdp_putprop(p,'sugar,s);
procedure vdp_setnumber(p,n);
vdp_putprop(p,'number,n);
procedure vdp_putprop(poly,prop,val);
begin scalar c,p;
c := cdr cdddr poly;
p := atsoc(prop,car c);
if p then
rplacd(p,val)
else
rplaca(c,(prop . val) . car c);
return poly
end;
procedure vdp_getprop(poly,prop);
(if p then cdr p else nil) where p=atsoc(prop,vdp_plist poly);
procedure vdp_fmon(coef,vev);
% Virtual distributive polynomial from monomial. [coef] is a BC;
% [vev] is EV. Returns a VDP, representing the monomial
% $[coef]^[vev]$
begin scalar r;
r := vdp_make(coef,vev,dip_fmon(coef,vev));
vdp_setsugar(r,ev_tdeg vev);
return r
end;
procedure vdp_2a(u);
dip_2a vdp_poly u;
procedure vdp_2f(u);
dip_2f vdp_poly u;
procedure vdp_2sq(u);
dip_2sq vdp_poly u;
procedure vdp_init(vars,sm,sx);
% Initializing vdp-dip polynomial package.
dip_init(vars,sm,sx);
procedure vdp_cleanup(l);
dip_cleanup(l);
procedure vdp_f2vdp(u);
begin scalar dip;
dip := dip_f2dip u;
if null dip then
return vdp_zero();
return vdp_make(dip_lbc dip,dip_evlmon dip,dip)
end;
procedure vdp_enumerate(f);
% f is a temporary result. Prepare it for medium range storage and
% assign a number.
if vdp_zero!? f or vdp_getprop(f,'number) then
f
else
vdp_putprop(f,'number,(vdp_pcount!* := vdp_pcount!* #+ 1));
procedure vdp_simpcont(p);
begin scalar q;
q := vdp_poly p;
if null q then
return p;
return vdp_fdip dip_simpcont q
end;
procedure vdp_lsort(pl);
% Distributive polynomial list sort. pl is a list of distributive
% polynomials. vdplsort(pl) returns the sorted distributive
% polynomial list of pl.
sort(pl,function vdp_evlcomp);
procedure vdp_evlcomp(p1,p2); % nicht auf das HM?
dip_evlcomp(vdp_poly p1,vdp_poly p2);
procedure vdp_ilcomb(v1,c1,t1,v2,c2,t2);
begin scalar r;
r := vdp_fdip dip_ilcomb(vdp_poly v1,c1,t1,vdp_poly v2,c2,t2);
vdp_setsugar(r,ev_max!#(
vdp_sugar v1 #+ ev_tdeg t1,vdp_sugar v2 #+ ev_tdeg t2));
return r
end;
procedure vdp_ilcombr(v1,c1,v2,c2,t2);
begin scalar r;
r := vdp_fdip dip_ilcombr(vdp_poly v1,c1,vdp_poly v2,c2,t2);
vdp_setsugar(r,ev_max!#(vdp_sugar v1,vdp_sugar v2 #+ ev_tdeg t2));
return r
end;
% The following procedures are not used for the Groebner basis
% computation but are useful anyway.
procedure gb_reduce(f,g);
% Groebner basis reduce. [f] is a VDP, [g] is a list of VDP's.
% Returns a VDP, a normal form of [f] wrt. [g]. Note that the
% result contains in general denominators.
begin scalar f1,c,vev,divisor,tai,fold;
fold := f;
f1 := vdp_setsugar(vdp_zero(),vdp_sugar f);
while not vdp_zero!? f do <<
vev := vdp_evlmon f;
c := vdp_lbc f;
divisor := gb_searchinlist(vev,g);
if divisor then <<
tai := T;
if vdp_monp divisor then
f := vdp_cancelmev(f,vdp_evlmon divisor)
else
f := gb_reduceonesteprat(f,c,vev,divisor);
>> else <<
f := gb_shift(f,f1);
f1 := secondvalue!*
>>
>>;
return if tai then f1 else fold
end;
procedure gb_reduceonesteprat(f,c,vev,g1);
% Groebner basis reduce one step rational. [f] is VDP; [c] is BC,
% the head coefficient of [f]; [vev] is a EV, the head term of [f];
% [g1] is a VDP such that its head monomial divides the head
% monomial of [f]. Returns a VDP $h$. $h$ is constructed from [f]
% by reducing the head monomial of [f] by [g1].
begin scalar b,rg1,vevcof;
rg1 := vdp_mred g1;
if vdp_zero!? rg1 then a % g1 is monomial
return vdp_mred f;
b := bc_quot(c,vdp_lbc g1);
vevcof := ev_dif(vev,vdp_evlmon g1);
return vdp_ilcombr(vdp_mred f,bc_a2bc 1,rg1,bc_neg b,vevcof);
end;
endmodule; % vdp
end; % of file