1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830 |
- module groebnr2; % Part 2 of the Groebner package.
-
- imports groebner,vdp2dip;
- create!-package('(groebnr2 groebman glexconv groebres groebmes
- groebrst groebtra groeweak hilberts hilbertp),
- '(contrib groebner));
- % Other packages needed.
- load!-package 'groebner;
- endmodule;
- module groebman; % Operators for manipulation of bases and
- % polynomials in Groebner style.
-
-
- fluid '(!*factor !*complex !*exp); % standard 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
- 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
- );
-
- 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);
-
- % variables for counting and numbering
- fluid '(hcount!* pcount!* mcount!* fcount!* bcount!* b4count!*
- basecount!* hzerocount!*);
-
- % control of the polynomial arithmetic actually loaded
- fluid '(currentvdpmodule!*);
- symbolic procedure gsorteval pars;
- % reformat a polynomial or a list of polynomials by a distributive
- % ordering; a list will be sorted and zeros are elimiated
- begin scalar vars,u,v,w,oldorder,nolist,!*factor,!*exp;
- integer n,pcount!*; !*exp := t;
- n := length pars;
- u := reval car pars;
- v := if n>1 then reval cadr pars else nil;
- if not eqcar(u,'list) then
- <<nolist := t; u := list('list,u)>>;
- w := for each j in groerevlist u
- collect if eqexpr j then !*eqn2a j else j;
- vars :=
- if null v then
- for each j in gvarlis w collect !*a2k j
- else groerevlist v;
- if not vars then vdperr 'gsort;
- oldorder := vdpinit vars;
- w := for each j in w collect a2vdp j;
- w := vdplsort w;
- w := for each x in w collect vdp2a x;
- while member(0,w) do w := delete(0,w);
- setkorder oldorder;
- return if nolist and w then car w else 'list . w;
- end;
-
- put('gsort,'psopfn,'gsorteval);
-
-
- symbolic procedure gspliteval pars;
- % split a polynomial into leading monomial and reductum;
- begin scalar vars,x,u,v,w,oldorder,!*factor,!*exp;
- integer n,pcount!*; !*exp := t;
- n := length pars;
- u := reval car pars;
- v := if n>1 then reval cadr pars else nil;
- u := list('list,u);
- w := for each j in groerevlist u
- collect if eqexpr j then !*eqn2a j else j;
- vars :=
- if null v then
- for each j in gvarlis w collect !*a2k j
- else groerevlist v;
- if not vars then vdperr 'gsplit;
- oldorder := vdpinit vars;
- w := a2vdp car w;
- if vdpzero!? w then x := w else
- <<x := vdpfmon(vdplbc w,vdpevlmon w);
- w := vdpred w>>;
- w := list('list,vdp2a x,vdp2a w);
- setkorder oldorder;
- return w;
- end;
- put('gsplit,'psopfn,'gspliteval);
- symbolic procedure gspolyeval pars;
- % calculate the S Polynomial from two given polynomials
- begin scalar vars,u,u1,u2,v,w,oldorder,!*factor,
- !*exp;
- integer n,pcount!*; !*exp := t;
- n := length pars;
- if n<2 or n#>3 then
- rerror(groebnr2,1,"GSpoly, illegal number or parameters");
- u1:= car pars; u2:= cadr pars;
- u := list('list,u1,u2);
- v := if n>2 then groerevlist caddr pars else nil;
- w := for each j in groerevlist u
- collect if eqexpr j then !*eqn2a j else j;
- vars := if null v then
- for each j in gvarlis w collect !*a2k j
- else v;
- if not vars then vdperr 'gspoly;
- oldorder := vdpinit vars;
- w := for each j in w collect a2vdp j;
- w := vdp2a groebspolynom3 (car w,cadr w);
- setkorder oldorder;
- return w;
- end;
-
- put('gspoly,'psopfn,'gspolyeval);
-
- symbolic procedure gvarseval u;
- % u is a list of polynomials; gvars extracts the variables from u
- begin integer n; scalar v,!*factor,!*exp; !*exp := t;
- n := length u;
- v := for each j in groerevlist reval car u collect
- if eqexpr j then !*eqn2a j else j;
- v := for each j in gvarlis v collect !*a2k j;
- v := if n= 2 then
- intersection (v,groerevlist reval cadr u) else v;
- return 'list . v
- end;
-
- put('gvars,'psopfn,'gvarseval);
-
-
- symbolic procedure greduceeval pars;
- % Polynomial reduction modulo a Groebner basis driver. u is an
- % expression and v a list of expressions. Greduce calculates the
- % polynomial u reduced wrt the list of expressions v reduced to a
- % groebner basis modulo using the optional caddr argument as the
- % order of variables.
- % 1 expression to be reduced
- % 2 polynomials or equations; base for reduction
- % 3 optional: list of variables
- begin scalar vars,x,u,v,w,np,oldorder,!*factor,!*groebfac,!*exp;
- integer n,pcount!*; !*exp := t;
- if !*groebprot then groebprotfile := list 'list;
- n := length pars;
- x := reval car pars;
- u := reval cadr pars;
- v := if n>2 then reval caddr pars else nil;
- w := for each j in groerevlist u
- collect if eqexpr j then !*eqn2a j else j;
- if null w then rerror(groebnr2,2,"Empty list in Greduce");
- vars :=
- if null v then
- for each j in gvarlis w collect !*a2k j
- else groerevlist v;
- if not vars then vdperr 'greduce;
- oldorder := vdpinit vars;
- groedomainmode();
- % cancel common denominators
- w := for each j in w collect reorder numr simp j;
- % optimize varable sequence if desired
- if !*groebopt then<< w:=vdpvordopt (w,vars); vars := cdr w;
- w := car w; vdpinit vars>>;
- w := for each j in w collect f2vdp j;
- if !*groebprot then w := for each j in w collect vdpenumerate j;
- if not !*vdpinteger then
- <<np := t;
- for each p in w do
- np := if np then vdpcoeffcientsfromdomain!? p
- else nil;
- if not np then <<!*vdpmodular:= nil; !*vdpinteger := t>>;
- >>;
- w := groebner2(w,nil);
- x := a2vdp x;
- if !*groebprot then
- <<w := for each j in w collect vdpenumerate j;
- groebprotsetq('candidate,vdp2a x);
- for each j in w do groebprotsetq(mkid('poly,vdpnumber j),
- vdp2a j)>>;
- w := car w;
- !*vdpinteger := nil;
- w := groebnormalform(x , w, 'sort);
- w := vdp2a w;
- setkorder oldorder;
- gvarslast := 'list . vars;
- return if w then w else 0;
- end;
-
-
- put('greduce,'psopfn,'greduceeval);
-
- % preduceeval moved to groesolv.red
- put('preduce,'psopfn,'preduceeval);
-
- endmodule;
- 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
- 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
- );
- 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,onlyheadtermreduction,groebprereduce,groebstat,
- gcheckpoint,gcdrart,gltbasis,groebsubs;
- % 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,
- 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,np,oldorder,!*gcd,!*ezgcd;
- 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;
- if not !*vdpinteger then
- np := t;
- for each p in w do
- np := np and vdpcoeffcientsfromdomain!? p;
- if not np then <<!*vdpmodular:= nil;
- !*vdpinteger := t;
- 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.vec)
- % such that cof * mon == vec 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,vec,r,f,fac1;
- if !*trgroeb then prin2t " ======= reducing ";
- fac1 := vdpgetprop(m,'factor);
- if fac1 then vec := vdpgetprop(fac1,'vector);
- if vec then
- << f := vdpprod(cdr vec,vdpgetprop(m,'monfac));
- cof := car vec>>
- 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;
- integer n;
- 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;
- module groebres;
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %
- % optimization of h-Polynomials by resultant calculation and
- % factorization
- %
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- fluid '(!*trgroeb );
- % the resultant is calculated from a h-polynomial and its predecessor
- % if both are bivariate in the same variables and if these variables
- % are the last ones in vdpvars*.
-
- expr procedure groebtestresultant(h1,h2,lv);
- begin scalar v1,hlist;
- v1 := indexcpl(vevsum0(lv,h1),1);
- if groebrescheck!?(2,v1,lv)
- and indexcpl(vevsum0(lv,h2),1)=v1
- then hlist :=
- reverse vdplsort
- groebhlistfromresultant
- (h1,h2,cadr reverse vdpvars!*)
- else
- if groebrescheck1!?(2,v1,lv)
- and indexcpl(vevsum0(lv,h2),1)=v1
- then hlist :=
- reverse vdplsort
- groebhlistfromresultant
- (h1,h2,caddr reverse vdpvars!*);
- if null hlist then return nil;
- return 'resultant .
- for each x in hlist collect list(h2,vdpenumerate x);
- end;
-
-
- expr procedure groebhlistfromresultant(h1,h0,x);
- % new h-polynomial calculation: calculate
- % the resultant of the two distributive polynomials h1 and h0
- % with respect to x.
- begin scalar ct00,hh,hh1,hs2;
- ct00:= time();
- hh:= vdpsimpcont groebresultant(h1,h0,x);
- if !*trgroeb then <<terpri();
- printb 57;
- prin2t " *** the resultant from ";
- vdpprint h1;
- prin2t " *** and";
- vdpprint h0;
- prin2t " *** is";
- vdpprint hh;
- printb 57;
- terprit 4 >>;
- hs2:= nil;
- if not vdpzero!? hh then
- << hh1:= vdp2a vdprectoint(hh,vdplcm hh);
- hh1:= factorf !*q2f simp hh1;
- if cdr hh1 and cddr hh1 then
- hs2:= for each p in cdr hh1 collect a2vdp prepf car p;
- if !*trgroeb and hs2 then
- <<prin2 " factorization of resultant successful:";
- terprit 2;
- for each x in hs2 do vdpprint x;
- terprit 2;
- ct00:= time() - ct00;
- prin2 " time for factorization:"; prin2 ct00;
- terpri() >>;
- >>;
- return hs2
- end;
- expr procedure groebresultant(p1,p2,x);
- begin scalar q1,q2,q;
- q1:=vdp2a vdprectoint(p1,vdplcm p1);
- q2:=vdp2a vdprectoint(p2,vdplcm p2);
- q:=a2vdp prepsq simpresultant list(q1,q2,x);
- return q;
- end;
- expr procedure groebrescheck!?(a,h1,vl);
- length h1 = a and car h1 = vl - 1;
- expr procedure groebrescheck1!?(a,h1,vl);
- length h1 = a and car h1 = vl - 2 and cadr h1 = vl - 1;
- endmodule;
- module groebmes;
-
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %
- % trace messages for the algorithms
- %
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- fluid '(vdpvars!* !*trgroeb !*trgroebs !*trgroeb1 !*trgroebr hcount!*
- hzerocount!* factorlevel!* basecount!* groetime!* pcount!*);
- symbolic procedure groebpairprint (p);
- <<groebmessff(" pair(",cadr p,nil);
- groebmessff(",",caddr p,nil);
- prin2 "), ";
- prin2 " lcm = ";print car p;
- >>;
- symbolic procedure groetimeprint;
- << prin2 ">> accum. cpu time:";
- prin2 (time() - groetime!*);
- prin2t " ms";
- >>;
- symbolic procedure groebmessff (m1,f,m2);
- << prin2 m1; prin2 vdpnumber f; if m2 then prin2t m2 >>;
- symbolic procedure groebmess0 (p);
- if !*trgroeb then << vdpprint p; >>;
- symbolic procedure groebmess1 (g,d);
- if !*trgroeb then <<
- prin2 " variables: "; print vdpvars!*;
- printbl();
- prin2t " Start of ITERATION "; terpri (); >>;
- symbolic procedure groebmess2 f;
- if !*trgroeb then << terpri();
- groebmessff ("polynomial ",f," eliminated");
- groetimeprint()>>;
- symbolic procedure groebmess2a(f,cf,fn);
- if !*trgroeb then << terpri();
- groebmessff ("polynomial ",f,nil);
- groebmessff (" elim. with cofactor ",cf," to");
- vdpprint fn; terpri();
- groetimeprint()>>;
- symbolic procedure groebmess3(p,s);
- if !*trgroebs then <<
- prin2 "S-polynomial from ";
- groebpairprint p;
- vdpprint s; terpri();
- groetimeprint();
- terprit 3 >>;
- symbolic procedure groebmess4 (p,d);
- << hcount!* := hcount!*+1;
- hzerocount!* := hzerocount!*+1;
- if !*trgroeb then <<
- terpri(); printbl();
- groebmessff(" reduction (",cadr p,nil);
- groebmessff(",", caddr p,nil);
- prin2 ") leads to 0; ";
- << prin2 n;
- prin2 if n=1 then " pair" else " pairs">>
- where n=length d;
- prin2t " left";
- printbl(); groetimeprint()>>;
- >>;
- symbolic procedure groebmess41 (p);
- << hcount!* := hcount!*+1;
- hzerocount!* := hzerocount!*+1;
- if !*trgroeb then << terpri(); printbl();
- groebmessff(" polynomial(", p,nil);
- prin2 ") reduced to 0;";
- terpri(); printbl(); groetimeprint()>>;
- >>;
-
- symbolic procedure groebmess5(p,h);
- if car p then % print for true h-Polys
- << hcount!* := hcount!* + 1;
- if !*trgroeb then << terpri();
- prin2 "H-polynomial ";
- prin2 pcount!*;
- prin2 " ev:"; prin2 vdpevlmon h;
- groebmessff(" from pair (",cadr p,nil);
- groebmessff(",", caddr p,")");
- vdpprint h; terpri();
- groetimeprint() >> >>
- else
- if !*trgroeb then << % print for input polys
- prin2t "from actual problem input:";
- vdpprint h;
- groetimeprint() >>;
- symbolic procedure groebmess50(g);
- if !*trgroeb1 then <<
- prin2 "list of active polynomials:";
- for each d1 in g do
- <<prin2 vdpgetprop(d1,'number);
- prin2 " ">>;
- terprit 2 >>;
- symbolic procedure groebmess51(d);
- if !*trgroeb1 then <<
- prin2t "Candidates for pairs in this step:";
- for each d1 in d do groebpairprint (d1);
- terprit 2 >>;
- symbolic procedure groebmess52(d);
- if !*trgroeb1 then <<
- prin2t "Actual new pairs from this step:";
- for each d1 in d do groebpairprint (d1);
- terprit 2 >>;
- symbolic procedure groebmess7 h;
- if !*trgroebs then
- <<prin2t "Testing factorization for "; vdpprint h>>;
- symbolic procedure groebmess8 (g,d);
- if !*trgroeb1 then <<
- prin2t " actual pairs: ";
- if null d then prin2t "null"
- else for each d1 in d do groebpairprint d1;
- groetimeprint() >>
- else if !*trgroeb then <<
- prin2 n; prin2t if n=1 then " pair" else " pairs" >>
- where n=length d;
- symbolic procedure groebmess13(g,problems);
- if !*trgroeb or !*trgroebr then <<
- if g then
- <<
- basecount!* := basecount!* +1;
- printbl(); printbl();
- prin2 "end of iteration ";
- for each f in reverse factorlevel!* do
- <<prin2 f; prin2 ".">>;
- prin2 "; basis "; prin2 basecount!*; prin2t ":";
- prin2 "{";
- for each g1 in g do vdpprin3t g1;
- prin2t "}";
- printbl(); printbl(); groetimeprint();
- >>
- else
- << printbl();
- prin2 "end of iteration branch ";
- for each f in reverse factorlevel!* do
- <<prin2 f; prin2 ".">>;
- prin2t " "; printbl(); groetimeprint();
- >>;
- if problems and !*trgroeb then
- <<
- groetimeprint(); terpri(); printbl();
- prin2 " number of partial problems still to be solved:";
- prin2t length problems;
- terpri();
- prin2 " preparing next problem ";
- if car car problems = 'file then
- prin2 cdr car problems
- else
- if cadddr car problems then
- vdpprint car cadddr car problems; % head of list G
- terpri();
- >> >>;
- symbolic procedure groebmess14 (h,hf);
- if !*trgroeb then <<
- prin2 "******************* factorization of polynomial ";
- (if x then prin2t x else terpri() )
- where x = vdpnumber(h);
- % vdpprint h;
- prin2t " factors:";
- for each g in hf do vdpprint car g;
- groetimeprint();
- >>;
- symbolic procedure groebmess15 f;
- if !*trgroeb then
- <<prin2t "***** monomial factor reduced:";
- vdpprint vdpfmon(a2vbc 1,f)>>;
- symbolic procedure groebmess19(p,restr,u);
- if !*trgroeb then <<
- printbl();
- prin2 "calculation branch ";
- for each f in reverse factorlevel!* do
- <<prin2 f; prin2 ".">>;
- prin2t " cancelled because";
- vdpprint p;
- prin2t "is member of an actual abort condition";
- printbl(); printbl();
- >>;
- symbolic procedure groebmess19a(p,u);
- if !*trgroeb then <<
- printbl();
- prin2 "during branch preparation ";
- for each f in reverse u do
- <<prin2 f; prin2 ".">>;
- prin2t " cancelled because";
- vdpprint p;
- prin2t "was found in the ideal branch";
- printbl();
- >>;
- symbolic procedure groebmess20 (p);
- if !*trgroeb then <<
- terpri();
- prin2 "secondary reduction starting with";
- vdpprint p>>;
- symbolic procedure groebmess21(p1,p2);
- if !*trgroeb then <<
- prin2 "polynomial ";
- prin2 vdpnumber p1;
- prin2 " replaced during secondary reduction by ";
- vdpprint p2;
- >>;
- symbolic procedure groebmess22(factl,abort1,abort2);
- if null factl then nil
- else
- if !*trgroeb then
- begin integer n;
- prin2t "BRANCHING after factorization point";
- n := 0;
- for each x in reverse factl do
- << n := n+1;
- prin2 "branch ";
- for each f in reverse factorlevel!* do
- <<prin2 f;prin2 ".";>>;
- prin2t n;
- for each y in car x do vdpprint y;
- prin2t "simple IGNORE restrictions for this branch:";
- for each y in abort1 do vdpprint y;
- for each y in cadr x do vdpprint y;
- if abort2 or caddr x then
- <<prin2t
- "set type IGNORE restrictions for this branch:";
- for each y in abort2 do vdpprintset y;
- for each y in caddr x do vdpprintset y;
- >>;
- printbl()>>;
- end;
- symbolic procedure groebmess23 (g0,rest1,rest2);
- if !*trgroeb then
- if null factorlevel!* then
- prin2t "** starting calculation ******************************"
- else
- << prin2 "** resuming calculation for branch ";
- for each f in reverse factorlevel!* do
- <<prin2 f; prin2 ".">>;
- terpri();
- if rest1 or rest2 then
- <<
- prin2t "-------IGNORE restrictions for this branch:";
- for each x in rest1 do vdpprint x;
- for each x in rest2 do vdpprintset x;
- >>;
- >>;
- symbolic procedure groebmess24(h,problems1,restr);
- % if !*trgroeb then
- <<prin2t
- "**********polynomial affected by NONNEGATIVE/POSITIVE condition:";
- vdpprint h;
- if restr then prin2t "under current restrictions";
- for each x in restr do vdpprint x;
- if null problems1 then prin2t " CANCELLED"
- else
- <<prin2t "partitioned into";
- vdpprintset car problems1;
- >>;
- >>;
- symbolic procedure groebmess25 (h,abort2);
- % if !*trgroeb then
- <<prin2t "reduction of set type cancel conditions by";
- vdpprint h;
- prin2t "remaining:";
- for each x in abort2 do vdpprintset x;
- >>;
- symbolic procedure groebmess26 (f1,f2);
- if !*trgroebs and not vdpequal(f1,f2) then
- <<terpri();
- prin2t "during final reduction";
- vdpprint f1;
- prin2t "reduced to";
- vdpprint f2;
- terpri();>>;
- symbolic procedure groebmess27 r;
- if !*trgroeb then
- <<terpri();
- prin2t "factor ignored (considered already):";
- vdpprint r>>;
- symbolic procedure groebmess27a (h,r);
- if !*trgroeb then
- <<terpri();
- vdpprint h;
- prin2t " reduced to zero by factor";
- vdpprint r>>;
- symbolic procedure groebmess28 r;
- if !*trgroeb then
- <<prin2 "interim content reduction:"; print r>>;
- symbolic procedure terprit n;
- % print blank lines.
- for i:=1:n do << terpri() >>;
- symbolic procedure printbl(); printb (linelength nil #- 2);
- symbolic procedure printb n; <<for i:=1:n do prin2 "-"; terpri()>>;
- symbolic procedure vdpprintset l;
- % print polynomials in one line;
- if l then
- << prin2 "{"; vdpprin2 car l;
- for each x in cdr l do
- <<prin2 "; "; vdpprin2 x>>;
- prin2t "}";>>;
- symbolic procedure vdpprin2l u;
- <<prin2 "("; vdpprin2 car u;
- for each x in cdr u do <<prin2 ","; vdpprin2 x;>>;
- prin2 ")";>>;
- endmodule;
- module groebrst;
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %
- % restrictions for polynomials during Groebner base calculation
- %
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- fluid '(groebrestriction!*);
-
- expr procedure groebtestrestriction (h,arg);
- if groebrestriction!* = 'nonnegative then groebnonneg(h,arg)
- else
- if groebrestriction!* = 'positive then groebpos(h,arg)
- else
- rerror(groebnr2,9,
- "Groebner: general restrictions not yet implemented");
- expr procedure groebnonneg(h,arg);
- % test if h is a polynomial which can have the value zero with
- % only nonnegative variable settings.
- begin scalar x,break,vev1,vevl,problems,problems1,r;
- if vdpzero!? h then return nil
- else
- if vevzero!? vdpevlmon h then goto finish;
- % car test the coefficients
- if vdpredzero!? h then return nil; % simple monomial
- break := nil;
- x := h;
- while not vdpzero!? x and not break do
- <<vev1 := vdpevlmon x;
- if not vbcplus!? vdplbc x then break := t;
- if not break then x := vdpred x>>;
- if break then return nil; % at least one negative coeff
- if vevzero!? vev1 then goto finish; % nonneg. solution imposs.
- % homogenous polynomial: find combinations of
- % variables which are solutions
- x := h;
- vev1 := vdpevlmon x;
- vevl := vevsplit(vev1);
- problems := for each x in vevl collect list x;
- x := vdpred x;
- while not vdpzero!? x do
- << vev1 := vdpevlmon x;
- vevl := vevsplit(vev1);
- problems1 := nil;
- for each e in vevl do
- for each p in problems do
- <<r := if not member(e,p) then e . p else p;
- problems1 := union(list r, problems1)>>;
- problems := problems1;
- x := vdpred x;
- >>;
- problems := % lift vevs to polynomials
- for each p in problems collect
- for each e in p collect
- vdpfmon(a2vbc 1,e);
- % rule out problems contained in others
- for each x in problems do
- for each y in problems do
- if not eq(x,y) and subset!?(x,y) then
- problems := delete (y,problems);
-
- % rule out some by cdr
- problems1 := nil;
- while problems do
- <<if vdpdisjoint!? (car problems,arg)
- then problems1 := car problems . problems1;
- problems := cdr problems;
- >>;
- finish:
- groebmess24(h,problems1,arg);
- return
- if null problems1 then 'cancel
- else 'restriction . problems1;
- end;
- expr procedure groebpos(h,dummy);
- % test if h is a polynomial which can have the value zero with
- % only positive (nonnegative and nonzero) variable settings.
- begin scalar x,break,vev1;
- if vdpzero!? h then return nil
- else
- if vevzero!? vdpevlmon h then return nil;
- % a simple monomial can never have pos. zeros
- if vdpredzero!? h then return groebposcancel(h);
- break := nil;
- x := h;
- % test coefficients
- while not vdpzero!? x and not break do
- <<vev1 := vdpevlmon x;
- if not vbcplus!? vdplbc x then break := t;
- if not break then x := vdpred x>>;
- if not break then return groebposcancel(h);
- if not groebposvevaluate h then groebposcancel(h);
- return nil;
- end;
- expr procedure groebposvevaluate h; t;
- % test if a polynomial can become zero under user restrictions
- % here a dummy to be rplaced elsewhere
- expr procedure groebposcancel(h);
- << groebmess24(h,nil,nil); 'cancel>>;
- endmodule;
- module groebtra;
-
- % calculation of a Groebner base with the Buchberger algorithm
- % including the backtracking information which denotes the
- % dependency between base and input polynomials
-
- % Authors: H. Melenk, H.M. Moeller, W. Neun
- % November 1988
-
- fluid '( % switches from the user interface
- !*groebopt !*groebfac !*groebres !*trgroeb !*trgroebs !*groebrm
- !*trgroeb1 !*trgroebr !*groebprereduce groebrestriction!*
- !*onlyheadtermreduction !*groebstat !*groebdivide !*groebprot
-
- vdpvars!* % external names of the variables
- !*vdpinteger !*vdpmodular % indicating type of algorithm
- 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.
- groetags!* % starting point of tag variables
- groetime!*
- );
-
- global '(gvarslast);
-
- switch groebopt,groebfac,groebrm,groebres,trgroeb,trgroebs,trgroeb1,
- trgroebr,onlyheadtermreduction,groebprereduce,groebstat,
- gcheckpoint,gcdrart,groebdivide,groebprot;
-
- % variables for counting and numbering
- fluid '(hcount!* pcount!* mcount!* fcount!* bcount!* b4count!*
- basecount!* hzerocount!*);
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % Interface
-
- symbolic procedure groebnertraeval u;
- % backtracking Groebner calculation
- begin integer n; scalar !*groebfac,!*groebrm,!*groebprot;
- n := length u;
- if n=1 then return groebnertra1(reval car u,nil,nil)
- else if n neq 2
- then rerror(groebnr2,10,
- "GROEBNERT called with wrong number of arguments")
- else return groebnertra1(reval car u,reval cadr u,nil)
- end;
- put('groebnert,'psopfn,'groebnertraeval);
- smacro procedure vdpnumber f;
- vdpgetprop(f,'number) ;
- symbolic procedure groebnertra1(u,v,mod);
- % Buchberger algorithm system driver. u is a list of expressions
- % and v a list of variables or NIL in which case the variables in u
- % are used.
- begin scalar vars,w,y,z,x,np,oldorder,groetags!*,tagvars;
- integer pcount!*,nmod;
- u := for each j in getrlist u
- collect
- <<if not eqcar(j,'equal) then 0 . j else
- if not idp cadr j then
- rerror(groebnr2,11,
- "GROEBNERT parameter 1 not in form {id = polynomial,...}")
- else cadr j . caddr j>>;
- if null u then rerror(groebnr2,12,"Empty list in Groebner")
- else if null cdr u then
- return 'list . list('equal,cdar u,caar u);;
- w := for each x in u collect cdr x;
- if mod then
- <<z := nil;
- for each vect in w do
- <<if not eqcar(vect,'list) then
- rerror(groebnr2,13,"Non list given to GROEBNER");
- if nmod=0 then nmod:= length cdr vect else
- if nmod<(x:=length cdr vect) then nmod=x;
- z := append(cdr vect,z);
- >>;
- if not eqcar(mod,'list) then
- rerror(groebnr2,14,"Illegal column weights specified");
- vars := groebnervars(z,v);
- tagvars := for i:=1:nmod collect mkid('! col,i);
- w := for each vect in w collect
- <<z:=tagvars; y := cdr mod;
- 'plus . for each p in cdr vect collect
- 'times . list('expt,car z,car y) .
- <<z := cdr z; y := cdr y;
- if null y then y := '(1); list p>>
- >>;
- groetags!* := length vars + 1;
- vars := append(vars,tagvars);
- >>
- else vars := groebnervars(w,v);
- groedomainmode();
- if null vars then vdperr 'groebner;
- oldorder := vdpinit vars;
- % cancel common denominators
- u := pair(for each x in u collect car x,w);
- u := for each x in u collect
- <<z := simp cdr x;
- multsq(simp car x,denr z ./ 1) . reorder numr z>>;
- w := for each j in u collect cdr j;
- % optimize varable sequence if desired
- if !*groebopt then<< w:=vdpvordopt (w,vars); vars := cdr w;
- w := car w; vdpinit vars>>;
- w := pair (for each x in u collect car x,w);
- w := for each j in w collect
- <<z := f2vdp cdr j; vdpputprop(z,'cofact,car j)>>;
- if not !*vdpinteger then
- <<np := t;
- for each p in w do
- np := if np then vdpcoeffcientsfromdomain!? p else nil;
- if not np then <<!*vdpmodular:= nil; !*vdpinteger := t>>;
- >>;
-
- w := groebtra2 w;
- w := if mod then groebnermodres(w,nmod,tagvars) else
- groebnertrares w;
- setkorder oldorder;
- gvarslast := 'list . vars;
- return w;
- end;
-
- symbolic procedure groebnertrares w;
- begin scalar c,u;
- return
- 'list . for each j in w collect
- <<c := prepsq vdpgetprop(j,'cofact);
- u := vdp2a j;
- if c=0 then u else list('equal,u,c)
- >>;
- end;
- symbolic procedure groebnermodres(w,n,tagvars);
- begin scalar x,c,oldorder;
- c := for each u in w collect prepsq vdpgetprop(u,'cofact);
- oldorder := setkorder tagvars;
- w := for each u in w collect
- 'list .
- <<u := numr simp vdp2a u;
- for i := 1:n collect
- prepf
- if not domainp u and mvar u = nth(tagvars,i)
- then <<x := lc u; u := red u; x>> else nil
- >>;
- setkorder oldorder;
- % reestablish term order for output
- w := for each u in w collect vdp2a a2vdp u;
- w := pair(w,c);
- return 'list .
- for each p in w collect
- if cdr p=0 then car p else
- list('equal,car p,cdr p);
- end;
-
-
-
- symbolic procedure preduceteval pars;
- % trace version of PREDUCE
- % parameters:
- % 1 expression to be reduced
- % formula or equation
- % 2 polynomials or equations; base for reduction
- % must be equations with atomic lhs
- % 3 optional: list of variables
- begin scalar vars,x,y,u,v,w,z,oldorder,!*factor,!*exp;
- integer n,pcount!*; !*exp := t;
- pars := groeparams(pars,2,3);
- y := car pars; u := cadr pars; v:= caddr pars;
- u := for each j in getrlist u
- collect
- <<if not eqcar(j,'equal) then
- j . j else cadr j . caddr j>>;
- if null u then rerror(groebnr2,15,"Empty list in PREDUCET");
- w := for each p in u collect cdr p; % the polynomials
- groedomainmode();
- vars :=
- if null v then
- for each j in gvarlis w collect !*a2k j
- else getrlist v;
- if not vars then vdperr 'preducet;
- oldorder := vdpinit vars;
- u := for each x in u collect
- <<z := simp cdr x;
- multsq(simp car x,denr z ./ 1) . reorder numr z>>;
-
- w := for each j in u collect
- <<z := f2vdp cdr j; vdpputprop(z,'cofact,car j)>>;
- if not eqcar(y,'equal) then y := list('equal,y,y);
- x := a2vdp caddr y; % the expression
- vdpputprop(x,'cofact,simp cadr y); % the lhs (name etc.)
- w := tranormalform(x,w, 'sort,'f);
- u := list('equal,vdp2a w,prepsq vdpgetprop(w,'cofact));
- setkorder oldorder;
- return u;
- end;
- put('preducet,'psopfn,'preduceteval);
-
- symbolic procedure groebnermodeval u;
- % Groebner for moduli calculation
- (if n=0 or n>3 then
- rerror(groebnr2,16,
- "GROEBNERM called with wrong number of arguments")
- else
- groebnertra1(reval car u,
- if n>=2 then reval cadr u else nil,
- if n>=3 then reval caddr u else '(list 1))
- ) where n = length u;
- put('groebnerm,'psopfn,'groebnermodeval);
-
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % some macros for local usage only
-
- smacro procedure tt(s1,s2);
- % lcm of leading terms of s1 and s2
- vevlcm(vdpevlmon s1,vdpevlmon s2);
-
-
-
- symbolic procedure groebtra2 (p);
- % setup all global variables for the Buchberger algorithm
- % printing of statistics
- begin scalar groetime!*,tim1,spac,spac1,p1,factortime!*,
- pairsdone!*,factorlvevel!*;
- factortime!* := 0;
- groetime!* := time();
- vdponepol(); % we construct dynamically
- hcount!* := pcount!* := mcount!* := fcount!* :=
- bcount!* := b4count!* := hzerocount!* := basecount!* := 0;
- if !*trgroeb then
- <<
- prin2 "Groebner Calculation with Backtracking starts ";
- terprit 2;
- >>;
- spac := gctime();
-
- p1:= groebtra3 (p);
-
- if !*trgroeb or !*trgroebr or !*groebstat then
- <<
- spac1 := gctime() - spac;
- terpri();
- prin2t "statistics for Groebner calculation";
- prin2t "===================================";
- prin2 " total computing time (including gc): ";
- prin2((tim1 := time()) - groetime!*);
- prin2t " milliseconds ";
- prin2 " (time spent for garbage collection: ";
- prin2 spac1;
- prin2t " milliseconds)";
- terprit 1;
- prin2 "H-polynomials total: "; prin2t hcount!*;
- prin2 "H-polynomials zero : "; prin2t hzerocount!*;
- prin2 "Crit M hits: "; prin2t mcount!*;
- prin2 "Crit F hits: "; prin2t fcount!*;
- prin2 "Crit B hits: "; prin2t bcount!*;
- prin2 "Crit B4 hits: "; prin2t b4count!*;
- >>;
- return p1;
- end;
-
-
- symbolic procedure groebtra3 (g0);
- begin scalar x,g,d,d1,d2,p,p1,s,h,g99,one;
- x := for each fj in g0 collect vdpenumerate trasimpcont fj;
- for each fj in x do g := vdplsortin(fj,g0);
- g0 := g; g := nil;
-
- % iteration :
- while (d or g0) and not one do
- begin
- if g0 then
- << % take next poly from input
- h := car g0; g0 := cdr g0;
- p := list(nil,h,h);
- >>
- else
- << % take next poly from pairs
- p := car d; d := cdr d;
- s := traspolynom (cadr p, caddr p); tramess3(p,s);
- h := groebnormalform(s,g99,'tree); %piloting wo cofact
- if vdpzero!? h then groebmess4(p,d)
- else h := trasimpcont tranormalform(s,g99,'tree,'h);
- >>;
- if vdpzero!? h then goto bott;
- if vevzero!? vdpevlmon h then % base 1 found
- << tramess5(p,h);
- g0 := d := nil; g:= list h;
- goto bott>>;
- s:= nil;
- % h polynomial is accepted now
- h := vdpenumerate h;
- tramess5(p,h);
- % construct new critical pairs
- d1 := nil;
- for each f in g do
- if groebmoducrit(f,h) then
- <<d1 := groebcplistsortin( list(tt(f,h),f,h) , d1);
- if tt(f,h) = vdpevlmon(f) then
- <<g := delete (f,g); groebmess2 f>>;
- >>;
- groebmess51(d1);
- d2 := nil;
- while d1 do
- <<d1 := groebinvokecritf d1;
- p1 := car d1;
- d1 := cdr d1;
- if groebbuchcrit4(cadr p1,caddr p1,car p1)
- then d2 := append (d2, list p1);
- d1 := groebinvokecritm (p1,d1);
- >>;
- d := groebinvokecritb (h,d);
- d := groebcplistmerge(d,d2);
- g := h . g;
- g99 := groebstreeadd(h, g99);
- groebmess8(g,d);
- bott:
- end; % ITERATION
- return groebtra3post g;
- end;
-
- symbolic procedure groebtra3post (g);
- % final reduction
- begin scalar r,p;
- g := vdplsort g;
- while g do
- <<p := tranormalform(car g,cdr g,'sort,'f);
- if not vdpzero!? p then r := trasimpcont p . r;
- g := cdr g>>;
- return reversip r;
- end;
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %
- % Reduction of polynomials
- %
-
- symbolic procedure tranormalform(f,g,type,mode);
- % general procedure for reduction of one polynomial from a set
- % f is a polynomial, G is a Set of polynomials either in
- % a search tree or in a sorted list
- % type describes the ordering of the set G:
- % 'TREE G is a search tree
- % 'SORT G is a sorted list
- % 'LIST G is a list, but not sorted
- % f has to be reduced modulo G
- % version for idealQuotient: doing side effect calculations for
- % the cofactors; only headterm reduction
- begin scalar c,vev,divisor,break;
- while not vdpzero!? f and not break do
- begin
- vev:=vdpevlmon f; c:=vdplbc f;
- divisor :=
- if type = 'tree then groebsearchinstree (vev,g)
- else groebsearchinlist (vev,g);
- if divisor and !*trgroebs then
- << prin2 "//-";
- prin2 vdpnumber divisor;
- >>;
- if divisor then
- if !*vdpinteger then
- f := trareduceonestepint(f,nil,c,vev,divisor)
- else
- f := trareduceonesteprat(f,nil,c,vev,divisor)
- else
- break := t;
- end;
- if mode = 'f then f := tranormalform1(f,g,type,mode);
- return f;
- end;
- symbolic procedure tranormalform1(f,g,type,mode);
- % reduction of subsequent terms
- begin scalar c,vev,divisor,break,f1;
- f1 := f;
- while not vdpzero!? f and not vdpzero!? f1 do
- <<f1 := f; break := nil;
- while not vdpzero!? f1 and not break do
- <<vev:=vdpevlmon f1; c:=vdplbc f1;
- f1 := vdpred f1;
- divisor :=
- if type = 'tree then groebsearchinstree (vev,g)
- else groebsearchinlist (vev,g);
- if divisor and !*trgroebs then
- << prin2 "//-";
- prin2 vdpnumber divisor;
- >>;
- if divisor then
- << if !*vdpinteger then
- f := trareduceonestepint(f,nil,c,vev,divisor)
- else
- f := trareduceonesteprat(f,nil,c,vev,divisor);
- break := t>>;
- >>;
- >>;
- return f;
- end;
-
-
-
- % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
- %
- % special reduction procedures
-
-
-
- symbolic procedure trareduceonestepint(f,dummy,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 vevlcm,a,b,cg,x,fcofa,gcofa;
- fcofa := vdpgetprop(f,'cofact);
- if null fcofa then fcofa := nil ./ 1;
- gcofa := vdpgetprop(g1,'cofact);
- if null gcofa then gcofa := nil ./ 1;
- vevlcm := vevdif(vev,vdpevlmon g1);
- cg := vdplbc g1;
- % calculate coefficient factors
- x := vbcgcd(c,cg);
- a := vbcquot(cg,x);
- b := vbcquot(c,x);
- f:= vdpilcomb1 (f, a, vevzero(),
- g1,vbcneg b, vevlcm);
- x := vdpilcomb1tra (fcofa, a, vevzero(),
- gcofa ,vbcneg b, vevlcm);
- vdpputprop(f,'cofact,x);
- return f;
- end;
-
-
- symbolic procedure trareduceonesteprat(f,dummy,c,vev,g1);
- % reduction step for rational case:
- % calculate f= f - g/vdpLbc(f)
- %
- begin scalar x,fcofa,gcofa,vev;
- fcofa := vdpgetprop(f,'cofact);
- gcofa := vdpgetprop(g1,'cofact);
- vev := vevdif(vev,vdpevlmon g1);
- x := vbcneg vbcquot(c,vdplbc g1);
- f := vdpilcomb1(f,a2vbc 1,vevzero(),
- g1,x,vev);
- x := vdpilcomb1tra (fcofa,a2vbc 1 , vevzero(),
- gcofa,x,vev);
- vdpputprop(f,'cofact,x);
- return f;
- end;
-
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %
- % calculation of an S-polynomial
- %
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- symbolic procedure traspolynom (p1,p2);
- begin scalar s,ep1,ep2,ep,rp1,rp2,db1,db2,x,cofac1,cofac2;
- if vdpzero!? p1 then return p1; if vdpzero!? p1 then return p2;
- cofac1 := vdpgetprop(p1,'cofact); cofac2 := vdpgetprop(p2,'cofact);
- ep1 := vdpevlmon p1; ep2 := vdpevlmon p2;
- ep := vevlcm(ep1, ep2);
- rp1 := vdpred p1; rp2 := vdpred p2;
- db1 := vdplbc p1; db2 := vdplbc p2;
- if !*vdpinteger then
- <<x:=vbcgcd(db1,db2); db1:=vbcquot(db1,x); db2:=vbcquot(db2,x)>>;
- ep1 := vevdif(ep,ep1); ep2 := vevdif(ep,ep2); db2 := vbcneg db2;
- s := vdpilcomb1 (rp2,db1,ep2,rp1,db2,ep1);
- x := vdpilcomb1tra (cofac2,db1,ep2,cofac1,db2,ep1);
- vdpputprop(s,'cofact,x);
- return s;
- end;
-
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %
- % Normalisation with cofactors taken into account
- %
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-
- symbolic procedure trasimpcont(p);
- if !*vdpinteger then trasimpconti p else trasimpcontr p;
-
- % routines for integer coefficient case:
- % calculation of contents and dividing all coefficients by it
-
- symbolic procedure trasimpconti (p);
- % calculate the contents of p and divide all coefficients by it
- begin scalar res,num,cofac;
- if vdpzero!? p then return p;
- cofac := vdpgetprop(p,'cofact);
- num := car vdpcontenti p;
- if not vbcplus!? num then num := vbcneg num;
- if not vbcplus!? vdplbc p then num:=vbcneg num;
- if vbcone!? num then return p;
- res := vdpreduceconti (p,num,nil);
- cofac:=vdpreducecontitra(cofac,num,nil);
- res := vdpputprop(res,'cofact,cofac);
- return res;
- end;
-
-
- % routines for rational coefficient case:
- % calculation of contents and dividing all coefficients by it
-
- symbolic procedure trasimpcontr (p);
- % calculate the contents of p and divide all coefficients by it
- begin scalar res,cofac;
- cofac := vdpgetprop(p,'cofact);
- if vdpzero!? p then return p;
- if vbcone!? vdplbc p then return p;
- res := vdpreduceconti (p,vdplbc p,nil);
- cofac:=vdpreducecontitra(cofac,vdplbc p,nil);
- res := vdpputprop(res,'cofact,cofac);
- return res;
- end;
-
- symbolic procedure vdpilcomb1tra (cofac1,db1,ep1,cofac2,db2,ep2);
- % the linear combination, here done for the cofactors
- % (standard quotients)
- addsq(multsq(cofac1,vdp2f vdpfmon(db1,ep1) ./ 1),
- multsq(cofac2,vdp2f vdpfmon(db2,ep2) ./ 1));
-
- symbolic procedure vdpreducecontitra(cofac,num,dummy);
- % divide the cofactor by a number
- quotsq(cofac,simp vbc2a num);
-
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %
- % special handling of moduli
- %
- symbolic procedure groebmoducrit(p1,p2);
- null groetags!*
- or pnth(vdpevlmon p1,groetags!*) = pnth(vdpevlmon p2,groetags!*);
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %
- % trace messages
- %
-
- symbolic procedure tramess0 x;
- if !*trgroeb then
- <<prin2t "adding member to intermediate quotient base:";
- vdpprint x;
- terpri()>>;
-
- symbolic procedure tramess1 (x,p1,p2);
- if !*trgroeb then
- <<prin2 "pair ("; prin2 vdpnumber p1; prin2 ",";
- prin2 vdpnumber p2;
- prin2t ") adding member to intermediate quotient base:";
- vdpprint x;
- terpri()>>;
-
-
- symbolic procedure tramess5(p,h);
- if car p then % print for true h-Polys
- << hcount!* := hcount!* + 1;
- if !*trgroeb then <<
- terpri();
- prin2 "H-polynomial ";
- prin2 pcount!*;
- groebmessff(" from pair (",cadr p,nil);
- groebmessff(",", caddr p,")");
- vdpprint h;
- prin2t "with cofactor";
- writepri(mkquote prepsq vdpgetprop(h,'cofact),'only);
- groetimeprint() >> >>
- else
- if !*trgroeb then << % print for input polys
- prin2t "candidate from input:";
- vdpprint h;
- groetimeprint() >>;
-
-
-
- symbolic procedure tramess3(p,s);
- if !*trgroebs then <<
- prin2 "S-polynomial ";
- prin2 " from ";
- groebpairprint p;
- vdpprint s;
- prin2t "with cofactor";
- writepri(mkquote prepsq vdpgetprop(s,'cofact),'only);
- groetimeprint();
- terprit 3 >>;
-
- endmodule;
- module groeweak; % weak test for f ~ 0 modulo G
- fluid '(!*groebweak current!-modulus pairsdone!* !*vdpinteger
- groebmodular!* !*groebfac);
- switch groebweak;
- symbolic procedure groebweakzerotest(f,g,type);
- % test f==0 modulo G with ON MODULAR
- begin scalar f1,c,vev,divisor,oldmode,a;
- integer n,oldprim;
- if vdpzero!? f then return f;
- if current!-modulus= 1 then setmod list 2097143;
- oldmode := setdmode('modular,t);
- f := groebvdp2mod f;
- f1 := vdpzero(); a:= vbcfi 1;
- while not vdpzero!? f and vdpzero!? f1 do
- begin
- vev:=vdpevlmon f; c:=vdplbc f;
- if type = 'sort then
- while g
- and vevcompless!? (vev,vdpevlmon (car g))
- do g := cdr g;
-
- divisor :=
- if type = 'tree then groebsearchinstree(vev,g)
- else groebsearchinlist (vev,g);
-
- if divisor and !*trgroebs then
- << prin2 "//m-";
- prin2 vdpnumber divisor;
- >>;
- if divisor then
- if vdplength divisor = 1 then
- f := vdpcancelmvev(f,vdpevlmon divisor)
- else
- <<divisor := groebvdp2mod(divisor);
- if divisor then f :=
- groebreduceonesteprat(f,nil,c,vev,divisor)
- else f1 := f>>
- else
- f1 := f;
- % ready:
- end;
- if not vdpzero!? f1 and !*trgroebs then
- <<prin2t " - nonzero result in modular reduction:";
- vdpprint f1;
- >>;
- setdmode('modular,nil);
- if oldmode then setdmode(get(oldmode,'dname),t);
- return vdpzero!? f1;
- end;
- smacro procedure tt(s1,s2);
- % lcm of leading terms of s1 and s2
- vevlcm(vdpevlmon s1,vdpevlmon s2);
- symbolic procedure groebweaktestbranch!=1(poly,g,d);
- % test GB(G) == {1} in modular style
- groebweakbasistest(list poly,g,d);
- symbolic procedure groebweakbasistest(g0,g,d);
- begin scalar oldmode,d,d1,d2,p,p1,s,h;
- scalar !*vdpinteger; % switch to field type calclulation
- return nil;
- if not !*groebfac then return nil;
- if current!-modulus= 1 then setmod list 2097143;
- if !*trgroeb then
- prin2t "---------------- modular test of branch ------";
- oldmode := setdmode('modular,t);
- g0 := for each p in g0 collect groebvdp2mod p;
- g := for each p in g collect groebvdp2mod p;
- d := for each p in d collect list (car p,
- groebvdp2mod cadr p, groebvdp2mod caddr p);
- while d or g0 do
- begin
- if g0 then
- << % take next poly from input
- h := car g0; g0 := cdr g0; p := list(nil,h,h);
- >>
- else
- << % take next poly from pairs
- p := car d;
- d := delete (p,d);
- s := groebspolynom (cadr p, caddr p);
- h:=groebsimpcontnormalform groebnormalform(s,g,'sort);
- if vdpzero!? h then !*trgroeb and groebmess4(p,d);
- >>;
-
- if vdpzero!? h then
- <<pairsdone!* := (vdpnumber cadr p . vdpnumber caddr p)
- . pairsdone!*;
- goto bott>>;
-
- if vevzero!? vdpevlmon h then % base 1 found
- << !*trgroeb and groebmess5(p,h);
- goto stop>>;
-
- s:= nil;
-
- h := vdpenumerate h; !*trgroeb and groebmess5(p,h);
- % construct new critical pairs
- d1 := nil;
- for each f in g do
- <<d1 := groebcplistsortin(list(tt(f,h),f,h),d1);
- if tt(f,h) = vdpevlmon(f) then
- <<g := delete (f,g);
- !*trgroeb and groebmess2 f>>;
- >>;
- !*trgroeb and groebmess51(d1);
- d2 := nil;
- while d1 do
- <<d1 := groebinvokecritf d1;
- p1 := car d1; d1 := cdr d1;
- d2 := groebinvokecritbuch4 (p1,d2);
- d1 := groebinvokecritm (p1,d1);
- >>;
-
- d := groebinvokecritb (h,d);
- d := groebcplistmerge(d,d2);
-
- g := h . g;
- goto bott;
-
- stop: d := g := g0 := nil;
-
- bott:
- end;
- if !*trgroeb and null g then
- prin2t "**** modular test detects empty branch!";
- if !*trgroeb then
- prin2t "------ end of modular test of branch ------";
- setdmode('modular,nil);
- if oldmode then setdmode(get(oldmode,'dname),t);
- return null g;
- end;
- fluid '(!*localtest);
- symbolic procedure groebfasttest(g0,g,d,g99);
- if !*localtest then
- <<!*localtest := nil;
- groebweakbasistest(g0,g,d)>>
- else if !*groebweak and g and vdpunivariate!? car g
- then groebweakbasistest(g0,g,d);
- symbolic procedure groebvdp2mod f;
- %convert a vdp in modular form
- % in case of headterm loss, nil is returned
- begin scalar u,c,mf;
- u := vdpgetprop(f,'modimage);
- if u then return if u='nasty then nil else u;
- mf := vdpresimp f;
- c := errorset!*(list('vbcinv,mkquote vdplbc mf),nil);
- if not pairp c then
- <<prin2t "************** nasty module (loss of headterm) ****";
- print f; print u; vdpprint f; vdpprint u;
- vdpputprop(f,'modimage,'nasty);
- return nil>>;
- u := vdpvbcprod(mf,car c);
- vdpputprop(u,'number,vdpgetprop(f,'number));
- vdpputprop(f,'modimage,u);
- return u;
- end;
- symbolic procedure groebmodeval(f,break);
- % evaluate LISP form r with REDUCE modular domain
- begin scalar oldmode,a,!*vdpinteger,groebmodular!*;
- groebmodular!* := t;
- if current!-modulus= 1 then setmod list 2097143;
- oldmode := setdmode('modular,t);
- a := errorset!*(f,t);
- setdmode('modular,nil);
- if oldmode then setdmode(get(oldmode,'dname),t);
- return if atom a then nil else car a;
- end;
- endmodule;
- module hilberts; % Hilbert series of a set of Monomials.
- % Author: Joachim Hollman, Royal Institute for Technology,
- % Stockholm, Sweden
- % email: <joachim@nada.kth.se>
- Comment
- --------------------------------------------------------
- A very brief "description" of the method used.
- M=k[x,y,z]/(x^2*y,x*z^2,y^2)
-
- x.
- 0 --> ker(x.) --> M --> M --> M/x --> 0
-
- M/x = k[x,y,z]/(x^2*y,x*z^2,y^2,x) = k[x,y,z]/(x,y^2)
-
- ker(x.) = ((x) + (x^2*y,x*z^2,y^2))/(x^2*y,x*z^2,y^2) =
-
- = (x,y^2)/(x^2*y,x*z^2,y^2)
-
-
- Hilb(ker(x.)) = Hilb - Hilb
- (x,y^2) (x^2*y,x*z^2,y^2)
-
- = 1/(1-t)^3 - Hilb -
- k[x,y,z]/(x,y^2)
-
- - (1/(1-t)^3 - Hilb
- k[x,y,z]/(x^2*y,x*z^2,y^2)
-
- = Hilb -Hilb
- M k[x,y,z]/(x,y^2)
-
- If you only keep the numerator in Hilb = N(t)/(1-t)^3
- M
-
-
- then you get
-
-
- (1-t)N (t) = N (t) - t (N (t) - N (t) )
- I I+(x) I Ann(x) + I
-
- i.e.
-
- N (t) = N (t) + t N (t) (*)
- I I+(x) Ann(x) + I
-
-
- Where
- I = (x^2*y,x*z^2,y^2)
- I + (x) = (x,y^2)
- I + Ann(x) = (x*y,z^2,y^2)
- N (t) is the numerator polynomial in Hilb
- I k[x,y,z]/I
-
-
- Equation (*) is what we use to compute the numerator polynomial, i.e.
- we "divide out" one variable at a time until we reach a base case.
- (One is not limited to single variables but I don't know of any good
- strategy for selecting a monomial.)
- Usage: hilb({monomial_1,...,monomial_n} [,variable])
- ;
- fluid '(nvars!*);
- % ************** MACROS ETC. **************
- smacro procedure term(c,v,e);
- list('times,c,list('expt,v,e));
- % -------------- safety check --------------
- smacro procedure var!-p(m);
- and(m,atom(m),not(numberp(m)));
- smacro procedure check!-expt(m);
- and(eqcar(m,'expt),var!-p(cadr(m)),numberp(caddr(m)));
- smacro procedure check!-single!-var(m);
- if var!-p(m) then t
- else check!-expt(m);
- smacro procedure check!-mon(m);
- if check!-single!-var(m) then t
- else if eqcar(m,'times) then check!-times(cdr(m))
- else nil;
- smacro procedure check!-args(monl,var);
- and(listp monl,eqcar(monl,'list),var!-p(var), check!-monl(monl));
- symbolic procedure make!-vector (n,pat);
- begin scalar v;
- v:=mkvect n;
- for i:=1:n do putv(v,i,pat);
- return v;
- end;
- % -------------- monomials --------------
- smacro procedure alloc!-mon(n);
- make!-vector(n,0);
-
- smacro procedure get!-nth!-exp(mon,n);
- getv(mon,n);
- smacro procedure set!-nth!-exp(mon,n,d);
- putv(mon,n,d);
- smacro procedure get!-tdeg(mon);
- getv(mon,0);
-
- smacro procedure set!-tdeg(mon,d);
- putv(mon,0,d);
- % -------------- ideals --------------
- smacro procedure the!-empty!-ideal();
- list(nil,nil);
- smacro procedure get!-next!-mon(ideal);
- <<x:=caadr ideal;
- if cdadr ideal then ideal:=list(car ideal,cdadr ideal)
- else ideal:=the!-empty!-ideal();
- x>>;
- smacro procedure not!-empty!-ideal(ideal);
- cadr ideal;
- smacro procedure first!-mon(ideal);
- caadr ideal;
- smacro procedure append!-ideals(ideal1,ideal2);
- list(car ideal2,append(cadr ideal1,cadr ideal2));
-
- symbolic procedure insert!-var(var,ideal);
- % inserts variable var as last generator of ideal
- begin
- scalar last;
- last:=list(make!-one!-var!-mon(var));
- return(list(last,append(cadr ideal,last)));
- end;
- symbolic procedure add!-to!-ideal(mon,ideal);
- % add mon as generator to the ideal
- begin
- scalar last;
- last := list(mon);
- if ideal = the!-empty!-ideal() then
- rplaca(cdr(ideal),last)
- else
- rplacd(car(ideal),last);
- rplaca(ideal,last);
- end;
- % ************** END OF MACROS ETC. **************
- % ************** INTERFACE TO ALGEBRAIC MODE **************
- symbolic procedure hilbsereval(u);
- begin
- scalar l,monl,var;
- l:=length(u);
- if l < 1 or l > 2 then
- rerror(groebnr2,17,
- "Usage: hilb({monomial_1,...,monomial_n} [,variable])")
- else if l = 1 then
- <<
- monl:= reval car u;
- var:= 'x;
- >>
- else
- <<
- monl:= reval car u;
- var:= reval cadr u;
- >>;
- monl := 'list . for each aa in (cdr monl) collect reval aa;
- if not check!-args(monl,var) then
- rerror(groebnr2,18,
- "Usage: hilb({monomial_1,...,monomial_n} [,variable])");
- % return(aeval
- % list('QUOTIENT,
- % coefflist2prefix(NPol(gltb2arrideal(monl)), var),
- % list('EXPT,list('PLUS,1,list('TIMES,-1,var)),
- % nvars!*)));
- return(aeval
- coefflist2prefix(npol(gltb2arrideal(monl)), var));
- end;
- % define "hilb" to be the algebraic mode function
- put('hilb,'psopfn,'hilbsereval);
- symbolic procedure check!-monl(monl);
- begin
- scalar flag,tmp;
- flag:=t;
- monl:=gltb!-fix(monl);
- while monl and flag do
- <<
- tmp:=car monl;
- flag:=check!-mon(tmp);
- monl:=cdr monl;
- >>;
- return(flag);
- end;
- symbolic procedure check!-times(m);
- begin
- scalar flag,tmp;
- flag:=t;
- while m and flag do
- <<
- tmp:= car m;
- flag:=check!-single!-var(tmp);
- m:=cdr m;
- >>;
- return(flag);
- end;
- symbolic procedure coefflist2prefix(cl,var);
- begin
- scalar i,poly;
- i:=0;
- for each c in cl do
- <<
- poly:= cons(term(c,var,i),poly);
- i:=i+1;
- >>;
- return (cons('plus,poly));
- end;
- symbolic procedure indets(l);
- % "indets" returns a list containing all the
- % indeterminates of l.
- % l is supposed to have a form similar to the variable
- % GLTB in the Groebner basis package.
- % (LIST (EXPT Z 2) (EXPT X 2) Y)
- begin
- scalar varlist;
- for each m in l do
- if m neq 'list then
- if atom(m) then varlist:=union(list(m),varlist)
- else if eqcar(m,'expt)
- then varlist:=union(list(cadr(m)),varlist)
- else varlist:=union(indets(cdr(m)),varlist);
- return(varlist);
- end;
- symbolic procedure build!-assoc(l);
- % Given a list of indeterminates (x1 x2 ...xn) we produce
- % an a-list of the form ((x1.1) (x2.2)... (xn.n))
- begin
- scalar i;
- i:=0;
- return(for each var in l collect progn(i:=i+1,cons(var,i)));
- end;
- symbolic procedure mons(l);
- % rewrite the leading monomials (i.e. GLTB).
- % the result is a list of monomials of the form:
- % (variable . exponent) or ((variable1 . exponent1) ...
- % (variablen . exponentn))
- %
- % mons('(LIST (EXPT Z 2) (EXPT X 2) (TIMES Y (EXPT X 3))));
- % (((Y . 1) (X . 3)) (X . 2) (Z . 2))
- begin
- scalar monlist;
- for each m in l do
- if m neq 'list then monlist:=
- if atom(m) then cons(cons(m,1),monlist)
- else if eqcar(m,'expt)
- then cons(cons(cadr(m),caddr(m)),monlist)
- else (for each x in cdr(m) collect mons!-aux(x)) . monlist;
- return(monlist);
- end;
- symbolic procedure mons!-aux(m);
- if atom(m) then cons(m,1)
- else if eqcar(m,'expt) then cons(cadr(m),caddr(m));
- symbolic procedure lmon2arrmon(m);
- % list-monomial to array-monomial
- % a list-monomial has the form: (variable_number . exponent)
- % or is a list with entries of this form.
- % "variable_number" is the number associated with the variable,
- % see build!-assoc()
- begin
- scalar tdeg,mon;
- mon:=alloc!-mon(nvars!*);
- tdeg:=0;
- if listp m then
- for each varnodotexp in m do
- <<
- set!-nth!-exp(mon,car varnodotexp, cdr varnodotexp);
- tdeg:=tdeg+cdr varnodotexp;
- >>
- else
- <<
- set!-nth!-exp(mon,car m, cdr m);
- tdeg:=tdeg+cdr m;
- >>;
- set!-tdeg(mon,tdeg);
- return(mon);
- end;
- symbolic procedure gltb!-fix(l);
- % sometimes GLTB has the form (list (list...))
- % instead of (list ...)
- if and(listp cadr l,caadr(l) = 'list) then cadr l else l;
-
- symbolic procedure ge(m1,m2);
- if get!-tdeg(m1) >= get!-tdeg(m2) then t else nil;
- symbolic procedure get!-end!-ptr(l);
- begin
- scalar ptr;
- while l do
- <<
- ptr:=l;
- l:=cdr l;
- >>;
- return(ptr);
- end;
- symbolic procedure gltb2arrideal(xgltb);
- % Convert the monomial ideal given by GLTB (in list form)
- % to a list of vectors where each vector represents a monomial.
- begin
- scalar l;
- l:=indets(gltb!-fix(xgltb));
- nvars!* := length(l);
- l:= sublis(build!-assoc(l),mons(gltb!-fix(xgltb)));
- l:=for each m in l collect lmon2arrmon(m);
- l:=sort(l,'ge);
- return(list(get!-end!-ptr(l),l));
- end;
-
- % ************** END OF INTERFACE TO ALGEBRAIC MODE **************
-
- %************** PROCEDURES **************
- symbolic procedure npol(ideal);
- % recursively computes the numerator of the Hilbert series
- begin
- scalar v,si;
- v:=next!-var(ideal);
- if not v then return(base!-case!-pol(ideal));
- si:=split!-ideal(ideal,v);
- return(shift!-add(npol(car si),npol(cadr si)));
- end;
-
- symbolic procedure divides!-by!-var(var,mon);
- begin
- scalar div;
- if get!-nth!-exp(mon,var) = 0 then return(nil);
- div:=alloc!-mon(nvars!*);
- for i:=1 : nvars!* do set!-nth!-exp(div,i,get!-nth!-exp(mon,i));
- set!-nth!-exp(div,var,get!-nth!-exp(mon,var)-1);
- set!-tdeg(div,get!-tdeg(mon)-1);
- return(div);
- end;
-
- symbolic procedure divides(m1,m2);
- % does m1 divide m2?
- % m1 and m2 are monomials
- % result: either nil (when m1 does not divide m2) or m2/m1
- begin
- scalar m,d,i;
- i:=1;
- m:=alloc!-mon(nvars!*);
- set!-tdeg(m,d:=get!-tdeg(m2)-get!-tdeg(m1));
- while d >= 0 and i <= nvars!* do
- <<
- set!-nth!-exp(m,i,d:=get!-nth!-exp(m2,i) - get!-nth!-exp(m1,i));
- i:= i+1;
- >>;
- if d < 0 then
- return (nil)
- else
- return(m);
- end;
- symbolic procedure shift!-add(p1,p2);
- % p1+z*p2
- % p1 and p2 are polynomials (nonempty coefficient lists)
- begin
- scalar p,pptr;
- pptr:=p:=cons(car p1,nil);
- p1:=cdr p1;
- while p1 and p2 do
- <<
- rplacd(pptr,cons(car p1 + car p2,nil));
- p1:=cdr p1;
- p2:=cdr p2;
- pptr:=cdr pptr;
- >>;
- if p1 then
- rplacd(pptr,p1)
- else
- rplacd(pptr,p2);
- return(p);
- end;
-
- symbolic procedure rem!-mult(ipp1,ipp2);
- % the union of two ideals with redundancy of generators eliminated
- begin
- scalar fmon,inew,isearch,primeflag,x;
- % fix; x is used in the macro...
- inew := the!-empty!-ideal();
- while not!-empty!-ideal(ipp1) and not!-empty!-ideal(ipp2) do
- begin
- if get!-tdeg(first!-mon(ipp2)) < get!-tdeg(first!-mon (ipp1))
- then <<
- fmon:=get!-next!-mon(ipp1);
- isearch:=ipp2;
- >>
- else
- <<
- fmon:=get!-next!-mon(ipp2);
- isearch:=ipp1;
- >>;
- primeflag:=t;
- while primeflag and not!-empty!-ideal(isearch) do
- if divides(get!-next!-mon(isearch),fmon) then primeflag:=nil;
- if primeflag then add!-to!-ideal(fmon,inew);
- end;
- if not!-empty!-ideal(ipp1) then return(append!-ideals(inew,ipp1))
- else return(append!-ideals(inew,ipp2));
- end;
- symbolic procedure next!-var(ideal);
- % extracts a variable in the ideal suitable for division
- begin
- scalar x,m,var;
- repeat
- << m:=get!-next!-mon(ideal);
- var:=get!-var!-if!-not!-single(m);
- >> until var or ideal = the!-empty!-ideal();
- return(var);
- end;
- symbolic procedure get!-var!-if!-not!-single(mon);
- % returns nil if the monomial is in a single variable,
- % otherwise the index of the second variable of the monomial
- begin
- scalar i,foundvarflag,exp;
- i := 0;
- foundvarflag:=nil;
- while not foundvarflag do
- <<
- i:=i+1;
- exp:=get!-nth!-exp(mon,i);
- if exp > 0 then foundvarflag:=t;
- >>;
- foundvarflag:=nil;
- while i < nvars!* and not foundvarflag do
- <<
- i:=i+1;
- exp:=get!-nth!-exp(mon,i);
- if exp > 0 then foundvarflag:=t;
- >>;
- if foundvarflag then return(i)
- else return(nil);
- end;
- symbolic procedure make!-one!-var!-mon(vindex);
- % returns the monomial consisting of the single variable vindex
- begin
- scalar mon;
- mon := alloc!-mon(nvars!*);
- for i := 1:nvars!* do set!-nth!-exp(mon,i,0);
- set!-nth!-exp(mon,vindex,1);
- set!-tdeg(mon,1);
- return(mon);
- end;
- symbolic procedure split!-ideal(ideal,var);
- % splits the ideal into two simpler ideals
- begin
- scalar div,ideal1,ideal2,m,x;
- ideal1 := the!-empty!-ideal();
- ideal2 := the!-empty!-ideal();
- while not!-empty!-ideal(ideal) do
- <<
- m:=get!-next!-mon(ideal);
- if div:=divides!-by!-var(var,m) then
- add!-to!-ideal(div,ideal2)
- else
- add!-to!-ideal(m,ideal1);
- >>;
- ideal2:=rem!-mult(ideal1,ideal2);
- ideal1:=insert!-var(var,ideal1);
- return(list(ideal1,ideal2));
- end;
- symbolic procedure base!-case!-pol(ideal);
- % in the base case every monomial is of the form Xi^ei
- % result : the numerator polynomial of the Hilbert series
- % i.e. (1-z^e1)*(1-z^e2)*...
- begin
- scalar p,degsofar,e,tdeg;
- tdeg:=0;
- for each mon in cadr ideal do tdeg:=tdeg + get!-tdeg(mon);
- p:=make!-vector(tdeg,0);
- putv(p,0,1);
- degsofar:=0;
- for each mon in cadr ideal do
- <<
- e:=get!-tdeg(mon);
- for j:= degsofar step -1 until 0 do
- putv(p,j+e,getv(p,j+e)-getv(p,j));
- degsofar:=degsofar+e;
- >>;
- return(vector2list(p));
- end;
- symbolic procedure vector2list v;
- % Convert a vector v to a list. No type checking is done.
- begin scalar u;
- for i:= upbv v step -1 until 0 do u := getv(v,i) . u;
- return u;
- end;
- endmodule;
- module hilbertp; % Computing Hilbert Polynomial from the Hilbert series.
- Comment
- Authors: H. Michael Moeller, Fernuniversitaet
- Hagen, Germany
- email: MA105@DHAFEU11.bitnet
- H. Melenk, Konrad-Zuse-Zentrum
- Berlin, Germany
- email: melenk@sc.ZIB-Berlin.de
- ;
- symbolic procedure newhilbi (bas,var,vars);
- begin scalar baslt,n,u,grad,h,joa,a,ii,dim0,varx;
- % extract leading terms
- baslt:= for each p in cdr bas collect
- <<u := hgspliteval list (p,vars); cadr u>>;
- % replace non atomic elements in the varlist by gensyms
-
- for each x in cdr vars do
- if (pairp x) then
- baslt := cdr subeval list(list('equal,x,gensym()),
- 'list . baslt);
- % compute the Hilbertseries
- joa := hilbsereval list ('list . baslt,var);
- % get the Hilbert polynomial
- grad := deg(joa,var);
- a:= for i:=0:grad collect coeffn (joa,var,i);
- n:= length cdr vars;
- %dim0 := (for i:=1:n product (var + i))/( for i:=1:n product i);
- varx := !*a2f var;
- dim0 := 1;
- for i:=1:n do dim0:= multf (addd(i,varx),dim0);
- dim0 := multsq(dim0 ./ 1,1 ./ (for i:=1:n product i));
- h := multsq(car(a) ./ 1,dim0);
- a := cdr(a);
- ii := 0;
-
- while a do
- << dim0 := multsq (dim0, addf(varx,numr simp (minus ii) )
- ./ addf(varx,numr simp (n -ii)));
- ii := ii + 1;
- if not (car a = 0) then h := addsq (h , multsq(car(a) ./ 1 ,dim0));
- a := cdr(a) >>;
- return mk!*sq h;
- end;
- symbolic procedure psnewhilbi (u);
- begin scalar zz;
- zz := 'list . gvarlis groerevlist reval car u;
- if length (u) = 2
- then return newhilbi ( reval car u, 'x, reval cadr u )
- else if length (u) = 1 then return newhilbi(reval car u,'x,zz)
- else rerror(groebnr2,19,"Wrong call to hilbertpolynomial");
- end;
-
- put ('hilbertpolynomial , 'psopfn ,'psnewhilbi);
- symbolic procedure hgspliteval pars;
- % A variant of Gsplit from grinterf.red.
- % Split a polynomial into leading monomial and reductum.
- begin scalar vars,x,u,v,w,oldorder,!*factor,!*exp;
- integer n,pcount!*; !*exp := t;
- n := length pars;
- u := reval car pars;
- v := if n>1 then reval cadr pars else nil;
- u := list('list,u);
- w := for each j in groerevlist u
- collect if eqexpr j then !*eqn2a j else j;
- vars :=
- if null v then
- for each j in gvarlis w collect !*a2k j
- else groerevlist v;
- if not vars then vdperr 'gsplit;
- oldorder := vdpinit vars;
- w := a2vdp car w;
- if vdpzero!? w then x := w else
- % <<x := vdpfmon(vdplbc w,vdpevlmon w);
- <<x := vdpfmon('( 1 . 1),vdpevlmon w);
- w := vdpred w>>;
- w := list('list,vdp2a x,vdp2a w);
- setkorder oldorder;
- return w;
- end;
-
- % simple Array access method for one- and two-dimensional arrays
- % NO check against misusage is done!
- % Usage: Rar:=makeRarray list dim1; Rar:=makeRarray list(dim1,dim2);
- % val:=getRarray(Rar,ind1); val:=getrarray(Rar,ind1,ind2);
- % putRarray(Rar,ind1,val); PutRarray(Rar,in1,ind2,val);
- % for two dimensional array access only !
- macro procedure functionindex2(u);
- begin scalar dims,ind1,ind2;
- dims := cadr u;
- ind1 := caddr u;
- ind2 := cadddr u;
- return %%%% ((ind1 #- 1) #* cadr dims) #+ ind2;
- list ('iplus,ind2,list('itimes,list('cadr,dims),
- list('iplus,ind1,-1)));
- end;
- macro procedure getrarray u;
- begin scalar arry,inds;
- arry := cadr(u);
- inds := cddr u;
- if length(inds) = 1 then
- return list('getv,list('cdr,arry),car inds)
- else
- return list('getv,list('cdr,arry),
- 'functionindex2 . list('car,arry) . inds);
- end;
- symbolic procedure makerarray dims;
- begin scalar u,n;
- n := for each i in dims product i;
- u := mkvect n;
- return dims . u;
- end;
- macro procedure putrarray u;
- begin scalar arry,inds, val;
- arry := cadr(u);
- inds := cddr u;
- val := nth (u,length u); % PSL: lastcar u;
- if length(inds) = 2 then
- return list('putv,list('cdr,arry),car inds,val)
- else
- return list('putv,list('cdr,arry), 'functionindex2 .
- list('car,arry).car inds. cadr inds . nil , val);
- end;
- procedure hilbertzerodimp(nrall, n, rarray);
- begin integer i, k, count, vicount;
- count := 0;
- i := 0;
- while ((i:= i+1) <= nrall and count < n) do
- begin vicount := 1;
- for k := 1:n do
- if (getrarray(rarray,i,k) = 0) then vicount := vicount + 1;
- if vicount = n then count := count + 1;
- end;
- return count = n;
- end;
-
- symbolic procedure groezerodim!?(f,n);
- begin scalar explist,a;
- integer r;
- %explist:= list( vev(lt(f1)),...,vev(lt(fr)) );
- explist:= for each fi in f collect vdpevlmon fi;
- r:= length(f);
- a := makerarray list(r,n);
- for i:=1 step 1 until r do
- for k:=1 step 1 until n do
- putrarray(a ,i,k, nth(nth(explist,i),k));
- return hilbertzerodimp (r, n, a);
- end;
- symbolic procedure gzerodimeval u;
- begin integer n;
- n := length u;
- if n = 2 then return gzerodim1(reval car u,reval cadr u)
- else
- rerror(groebnr2,20,
- "GZERODIM? called with wrong number of arguments")
- end;
- put('gzerodim!?,'psopfn,'gzerodimeval);
- symbolic procedure gzerodim1(u,v);
- begin scalar vars,w,z,dv,oldorder;
- w := for each j in getrlist u
- collect if eqexpr j then !*eqn2a j else j;
- if null w then rerror(groebnr2,21,
- "Empty list in HilbertPolynomial");
- vars :=
- if null v then for each j in gvarlis(w) collect !*a2k j
- else % test, if vars are really used
- << z := gvarlis (w);
- for each j in (v:= getrlist v) do
- if member(j,z) then dv := !*a2k j . dv;
- dv := reversip dv;
- if not (length v = length dv) then
- << prin2 " HilbertPolynomial: ";
- prin2 (length v - length dv);
- prin2t " of the variables not used";
- terpri () >>;
- dv>>;
- if not vars then vdperr 'gzerodim!?;
- oldorder := vdpinit vars;
- w := for each j in w collect numr simp j;
- w := for each j in w collect f2vdp j;
- w := groezerodim!? (w,length vars);
- setkorder oldorder;
- return if w then newhilbi(u,'x,'list . v) else nil;
- end;
- symbolic procedure gbtest(g);
- % test, if the given set of polynomials is a Groebner basis.
- % only fast to compute plausilbility test.
- begin scalar fredu,g1,r,s;
- g := vdplsort g;
- % make abbreviated version of G
- g1:= for each p in g collect
- << r := vdpred p;
- if vdpzero!? r then p else
- vdpsum(vdpfmon(vdplbc p,vdpevlmon p),
- vdpfmon(vdplbc r,vdpevlmon r))
- >>;
- while g1 do
- <<for each p in cdr g1 do
- if not groebbuchcrit4t(vdpevlmon car g1,vdpevlmon p) then
- <<s := groebspolynom(car g1,p);
- if not vdpzero!? s and
- null groebsearchinlist (vdpevlmon s,cddr g1)
- then rerror(groebnr2,22,
- "****** Not a GROEBNER basis wrt current ordering");
- >>;
- if groebsearchinlist(vdpevlmon car g1,cdr g1) then
- fredu := t;
- g1 := cdr g1;
- >>;
- if fredu then
- <<terpri!* t;
- prin2t "WARNING:system is not a fully reduced GROEBNER basis";
- prin2t "with current term ordering";
- >>;
- end;
-
- endmodule;
- end;
|