module groesolv; % Tools for solving systems of polynomials (and poly-
% % nomial equations) based on Groebner basis techniques.
% Authors: H.Melenk (ZIB Berlin)
% H.M Moeller (Fernunversitaet Hagen)
% W.Neun (ZIB Berlin)
%
% Aug 1992: accepting external solutions for univariate pols.
%
% March 1994: access to roots/multroot.
%
% Feb 1995: assumptions, requirements added.
% operators:
%
% GROESOLVE does the whole job of solving a nonlinear
% set of expressions and / or equations
%
% GROEPOSTPROC expects that its car parameter is a
% lexical groebner base already
fluid '(groesolvelvevel!* !*groebprereduce);
fluid '(groesoldb!* groesoldmode!* !*groebnumval !*groebcomplex
!*groebopt !*groebprot !*groesolrecurs
!*groesolgarbage denominators!* variables!*
groebroots!* % a-list for predefined root expressions
depl!* % the reduce dependency list
!*vdpinteger % coefficient mode
!*compxroots % switch for multroot
gmodule
!*arbvars
);
fluid '(!*convert !*ezgcd !*msg !*precise dmode!* !*varopt);
global '(groebprotfile gvarslast !*trgroesolv glterms assumptions
requirements);
groesolvelvevel!* := 0;
symbolic procedure groesolveeval u;
begin scalar !*ezgcd,GBlist,oldtorder,res,!*groebopt,!*precise,
y,fail,problems,denominators!*,variables!*,gmodule,at;
if null dmode!* then !*ezgcd := t;
!*groebopt := !*varopt; gvarslast :='(list);
oldtorder := apply1('torder,'(lex));
groesoldmode!* := get(dmode!*,'dname);
!*groebcomplex := !*complex;
groesetdmode(groesoldmode!*,NIL);
problems:={u};
while problems and not fail do
<<u:=car problems; problems:=cdr problems;
GBlist := cdr GroebnerFeval u;
at := union(at,cdr glterms);
!*groebopt := nil; % 29.8.88
groesetdmode(groesoldmode!*,T);
if null variables!* then variables!*:=gvarslast;
if not(Gblist = '((list 1))) then
for each gb in GBlist do
<<if !*trgroesolv then
<<writepri("starting from basis",'first);
writepri(mkquote gb,'last)>>;
for each r in res do if gb
% do not compare with the mother problem;
and not subsetp(car r,car u)
then
if groesolvidsubset!?(gb,car r,variables!*) then
res:=delete(r,res) else
if groesolvidsubset!?(car r,gb,variables!*) then
<<gb:=nil;
if !*trgroesolv then writepri("redundant",'only)
>>;
if gb then
<<y:=groesolvearb(groesolve0(gb,variables!*),variables!*);
if y neq 'failed then res:=(gb.y).res else fail:=t;
if !*trgroesolv then
<<writepri("partial result:",'first);
writepri(mkquote('list.cdar res),'last)
>>;
for each d in denominators!* do
problems:={append(gb,{d}),variables!*}.problems;
denominators!*:=nil;
>>;
>>;
>>;
apply1('torder,{oldtorder});
problems:=nil; if fail then res:=nil;
if null res then requirements:=append(requirements,at)
else assumptions := append(assumptions,at);
for each r in res do problems:=union(cdr r,problems);
return 'list . groesolve!-redun2 problems
end;
symbolic procedure groesolve!-redun2 sol;
% Sol is a list of solutions; remove redundancies, now not
% by ideal theory but by simple substitution.
begin scalar b;
for each s in sol do
if s memq sol then
<<b:=nil;
for each r in sol do
if not(r eq s) then
if not b and groesolve!-redun2a(r,s) then b:=r;
if b then sol:=delete(s,sol);
>>;
return sol;
end;
symbolic procedure groesolve!-redun2a(r,s);
% redundancy test: if sub(s,r)=> trivial then t because s
% is a special case of r.
if smemq('root_of,s) then nil else
begin scalar q,!*evallhseqp,!*protfg;
!*evallhseqp:=t; !*protfg:=t;
q:=errorset({'subeval,mkquote {s,r}},nil,nil);
if errorp q then <<erfg!*:=nil;return nil>>;
q:=cdar q;
while q and 0=reval{'difference,cadar q,caddar q} do q:=cdr q;
return null q;
end;
symbolic procedure groesolvidsubset!?(b1,b2,vars);
% test if ideal(b1) is a subset of ideal(b2)
% (b2 is a specialization of b1 wrt zeros).
null b1 or (car b1='list or 0=preduceeval{car b1,b2,vars}) and
groesolvidsubset!?(cdr b1,b2,vars);
symbolic procedure groesolvearb(r,vars);
% Cover unmentioned variables.
if atom r or not !*arbvars then r else
for each s in r collect
<<for each x in cdr vars do
if not smember(x,s) then
s:=append(s,{{'equal,x,prepf makearbcomplex()}});
s>>;
%------------------- driver for the postprocessor ----------------
symbolic procedure groesolve0(A,vars);
begin scalar r,ids,newvars,newA;
if(r:=groepostnumsolve(A,vars))then return r;
if(r:=groepostfastsolve(A,vars))then return r;
r:=groepostsolveeval{A,vars};
if r neq 'failed then return cdr r;
ids:=cdr gindependent_seteval{A,vars};
if null ids then goto nullerr;
ids:=car ids;
newvars:='list.for each x in cdr vars join
if not(x memq ids) then {x};
newA:=groebnereval{A,newvars};
denominators!*:=cdr glterms;
if newA='(list 1) then rerror(groebner,24,
"recomputation for dim=0 failed");
r:=groepostfastsolve(newA,newvars);
if r then return r;
r:=groepostsolveeval{A,vars};
if r neq 'failed then return cdr r;
nullerr:
rerror(groebner,23,
"Moeller ideal decomposition failed with 0 dim ideal");
end;
symbolic procedure groepostnumsolve(gb,vars);
if not errorp errorset('(load!-package 'roots2),nil,nil)
and getd 'multroot0
and get(dmode!*,'dname) member '(ROUNDED COMPLEX!-ROUNDED)
and length gb = length vars and groepostnumsolve1(gb,vars)
then (cdr reval multroot0(precision 0,gb)) where !*compxroots=t;
symbolic procedure groepostnumsolve1(gb,vars);
if null gb then t else
groepostnumsolve1(cdr gb,cdr vars) and
<<for each x in kernels numr simp car gb do q:=q and x member vars;
q>> where q=t;
symbolic procedure groepostfastsolve(gb,vars);
% try to find a fast solution.
begin scalar u,p1,p2,fail,kl,res;
if !*trgroesolv then prin2t "fast solve attempt";
groesoldmode!* := get(dmode!*,'dname);
!*groebnumval := member(groesoldmode!* ,
'(ROUNDED COMPLEX!-ROUNDED));
groesetdmode(groesoldmode!*,'NIL);
u:=kl:=for each p in cdr gb collect
<<p:=numr simp p;
intersection(vars,kernels p).p>>;
if u='((nil)) then goto trivial;
while u and cdr u do
<<p1:=car u; p2:=cadr u; u:= cdr u;
car p1:=setdiff(car p1,car p2);
fail:=fail or null car p1>>;
if fail then goto exit;
res:= for each r in groepostfastsolve1(reverse kl,nil,0)
collect 'list.reverse r;
goto exit;
trivial:
res:= {'list.for each x in cdr vars collect
{'equal,x,mvar makearbcomplex()}};
exit:
groesetdmode(groesoldmode!*,t);
return res;
end;
fluid '(f);
symbolic procedure groepostfastsolve1(fl,sub,n);
if null fl then '(nil) else
begin scalar u,f,v,sub1;
n:=n+1;
f:=car fl; v:=car f; f:=numr subf(cdr f,sub);
if null f then return groepostfastsolve1(cdr fl,sub,n);
% v:=car sort(v,function(lambda(x,y);degr(f,x)>degr(f,y)));
v := car v;
(f:=reorder f) where kord!*={v};
if not domainp lc f then groepostcollectden reorder lc f;
u:=groesolvePolyv(prepf f,v);
return for each s in u join
<<sub1:=if smemq('root_of,s) then sub else
(v . caddar s) . sub;
for each q in groepostfastsolve1(cdr fl,sub1,n) collect
car s.q
>>;
end;
unfluid '(f);
symbolic procedure groepostcollectden d;
% d is a non trivial denominator (standard form);
% collect its factors.
for each p in cdr fctrf d do
if not member(p:=prepf car p,denominators!*) then
denominators!*:=p.denominators!*;
put('groesolve,'psopfn,'groesolveeval);
symbolic procedure groepostsolveeval u;
begin scalar a,b,vars,oldorder;
scalar groesoldb!*;
scalar !*groebprereduce,!*groebopt,!*groesolgarbage;
groesoldmode!* := get(dmode!*,'dname);
groesetdmode(groesoldmode!*,'NIL);
!*groebnumval := member(groesoldmode!* ,
'(ROUNDED COMPLEX!-ROUNDED));
if vdpsortmode!* = 'LEX then t else
rerror(groebner,8,
"Groepostproc, illegal torder; (only LEX allowed)");
a := groerevlist reval car u;
vars := cdr u and groerevlist cadr u or groebnervars(a,nil);
oldorder := setkorder vars;
b := groesolve1(a,a,vars);
a := NIL;
if b eq 'failed then a:=b else
<<for each sol in b do % delete multiplicities
if not member(sol,a) then a := sol . a;
a := 'list . for each sol in a collect 'list . sol;
>>;
setkorder oldorder;
groesetdmode(groesoldmode!*,t);
return a;
end;
put('groepostproc ,'psopfn,'groepostsolveeval);
% DATA structure:
%
% all polynomials are held in prefix form (expressions)
% transformation to standard quotients/ standard forms is done locally
% only; distributive form is not used here.
%
% a zero is a set of equations, if possible with a variable on the
% lhs each
% e.g. {y=17,z=b+8};
% internally: ((equal y 17)(equal z (plus b 8)))
% a zeroset is a list of zeros
% elgl {{y=17,z=b+8},{y=17,z=b-8}}
% internally the sets (lists) are kept untagged as lists; the
% tag 'list is only added to the results and to those lists which
% are parameters to algebraic operators not in this package.
symbolic procedure groesolve1 (A,fullA,vars);
% A lex Groebner basis or tail of lex Groebner basis
% fullA the complete lex Groebner basis to A
% vars the list of variables
if null A or A='(1) then nil else
<<begin scalar f1,A1,res,Q,gi,Ng1,Ng2,NgQ,Qg,mv,test;
res := assoc (A,groesoldb!*); if res then return cdr res;
groesolvelvevel!* := groesolvelvevel!* + 1;
%%tr prin2t "=================================================";
%%tr prin2l list( " groesolvelvevel!*: ",groesolvelvevel!*);
%%tr prin2t " Problem:";
%%tr writepri (mkquote ('list . a),'only);;
if member(A,!*groesolrecurs) then return 'failed;
% <<vars := append(cdr vars,{car vars});
% if member(vars,!*groesolrecurs) then
% <<!*groesolgarbage := T;
% return list for each p in A collect list('EQUAL,p,0) >>;
% !*groesolrecurs := vars . !*groesolrecurs;
% A := cdr Groebnereval{'list.A,'list . vars};
% >>;
!*groesolrecurs := A . !*groesolrecurs;
if length A = 1 then
<< %%tr print "single polynomial; solve it";
res := groesolvePoly car A; goto ready>>;
% step 1
f1 := car A;
A1 := cdr A;
test := nil;
mv := intersection(vars,LtermVariables f1); % test Buchcrit 4
for each p in A1 do
if intersection (mv,LtermVariables p) then test := T;
if not test then
<< %%tr print "f1 has unique main variable";
NgQ := groesolve1(A1,A1,vars);
if NgQ eq 'failed then <<res:='failed;goto ready>>;
res := ZerosetIntersection(NgQ,f1,vars);
goto ready>>;
% Q := cdr GroebIdq('list . A1,f1,'list . vars); % A1:f1
Q := GroesolvIdq(A1,f1,vars);
%%tr print "A1:f1"; %%tr writepri (mkquote ('list . Q),'only);;
if Q='(1) then % f1 already member of A1; skip it
<<%%tr print "f1 already in A1; ignore";
res:= groesolve1(A1,fullA,vars);
goto ready>>;
NgQ := groesolve1(Q,Q,vars);
if NgQ eq 'failed then <<res:='failed;goto ready>>;
ng1 := ZerosetIntersection (NgQ,f1,vars);
% step 4
if GroesolvIdqTest(A1,f1,vars) then
<< while Q do
<< gi := car Q; Q := cdr Q;
gi := preduceeval list (gi , 'list . A, 'list . vars);
if gi = 0 then Q:= nil else
A := cdr GroebIdq('list . A ,gi,'list . vars);
>>;
Ng2 := groesolve1(A,A,vars);
if Ng2 eq 'failed then <<res:='failed;goto ready>>;
>> else
<<Ng2 := ();
if length Q = 1 then
<< gi := preduceeval list (car Q, 'list . fullA, 'list . vars);
if gi neq 0 then
<< Qg := cdr GroebIdq('list . fullA,gi,'list . vars); % A1:gi
Ng2 := groesolve1(Qg,Qg,vars);
if Ng2 eq 'failed then <<res:='failed;goto ready>>;
>> >>
else
<<Ng2 := groesolve2(A1,Q,vars);
if Ng2 eq 'failed then <<res:='failed;goto ready>>
>>
>>;
res:= ZerosetUnion (Ng1,Ng2);
ready:
%%tr print list( "Result, level!*: ",groesolvelvevel!*);
%%tr writeres res;
%%tr print "...................................................";
groesolvelvevel!* := groesolvelvevel!* -1;
groesoldb!* := (a . res) . groesoldb!*;
return res;
end
>> where !*groesolrecurs = !*groesolrecurs; % recursive fluid!
symbolic procedure GroesolvIdqTest(A1,f1,vars);
not(deg(f1,car vars) eq deg(car A1,car vars));
symbolic procedure GroesolvIdq(A1,f1,vars);
begin scalar x,temp;
x := car vars;
if not GroesolvIdqTest(A1,f1,vars)
then return cdr GroebIdq('list . A1,f1,'list . vars);
temp :=
for each p in A1 collect
reval car reverse coeffeval list(p,x);
return cdr groebnereval list ('list . temp,'list . vars);
end;
symbolic procedure groesolve2(A,Q,vars);
% Calculation of the zeroset A1:(g1,g2,,,,gs),
% the gi given as members of Q.
groesolvetree (A,Q,Q,vars);
symbolic procedure groesolvetree(A,T1,phi,vars);
begin scalar Q,NGtemp,NGall,T2,h,g,A2,phi2;
%%tr prin2t ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>";
%%tr prin2t "groesolvetree; A:";
%%tr writepri(mkquote ('list . a),'only);
%%tr writepri( "T1 :",'car);
%%tr writepri(mkquote ('list . T1) ,'last);
%%tr writepri( "phi:",'car);
%%tr writepri(mkquote ('list . phi),'last);
if null phi then return nil else
if null T1 then
<<Q := cdr GroebIdq('list.A, 'TIMES.phi,'list.vars);
%%tr prin2t "<< << << << << << << << << << << << << <<<";
return if car Q = 1 then nil
else groesolve1(Q,Q,vars) >>;
for each g in T1 do
<<h:=Preduceeval list(g,'LIST.A,'LIST.vars);
phi := delete(g,phi);
if not(h=0) then <<T2:=h.T2; phi:=h.phi>>;
>>;
if null phi then return nil; % 29.8.88
T1 := T2;
Q := cdr GroebIdq('list.A, 'TIMES.phi,'list.vars);
%%tr writepri( "PHI :",'car);
%%tr writepri(mkquote ('TIMES.phi) ,'last);
%%tr writepri( "A:PHI :",'car);
%%tr writepri(mkquote ('LIST.Q) ,'last);
if not(car Q = 1) then
<< NGall := groesolve1(Q,Q,vars);
if NGall eq 'failed then return 'failed;
>>;
if !*groesolgarbage then
return groesolverestruct(Q,phi,vars,NGall);
while T1 do
<<g:=car T1; T1:=cdr T1;
phi2 := delete(g,phi);
if phi2 then
<<A2 := cdr Groebnereval list('LIST . g . A,'LIST . vars);
if not(car A2 =1) then
<<NGtemp := groesolvetree(A2,T1,phi2,vars);
NGall := ZerosetUnion(NGtemp,NGall)>>;
>>>>;
%%tr prin2t "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<";
return NGall;
end;
symbolic procedure groesolverestruct(A,phi,vars,NGall);
% there was a problem with an embedded solution in phi such that
% A:phi = A
% we try a heuristic by making one variable a formal parameter
begin scalar newA,newvars,mv,oldorder,solu;
mv := LtermVariables ('TIMES.phi);
mv := car mv;
%%tr prin2 "groesolverestruct:taking variable ";prin2t mv;
newvars := delete(mv,vars);
oldorder := setkorder newvars;
newA := cdr GroebnerEval list('LIST . A,'LIST . newvars);
%%tr prin2t "new Groebner basis:";
%%tr writepri(mkquote ('LIST . newA),'ONLY);
!*groesolgarbage := NIL;
solu := groesolve1(newA,newA,newvars);
%%tr if !*groesolgarbage then prin2t "Heuristics failed"
%%tr else prin2t "better result";
setkorder oldorder;
return if !*groesolgarbage then NGall else solu;
end;
%%tr symbolic procedure writeres r;
%%tr writepri (mkquote ('list.for each x in r collect 'list.x),'ONLY);
symbolic procedure LtermVariables u;
% extract variables of leading term in u
begin scalar v;
u := numr simp u;
while not domainp u do
<<v := mvar u . v;
u := lc u>>;
return reversip v;
end;
symbolic procedure ZerosetIntersection (NG,poly,vars);
% NG is a zeroset
% poly is a polynomial
% the routine maps the zeros in NG by the polynomial:
% each zero is substituted into the polynomial;
% that gives a univariate
% solved by SOLVE or numerical techniques.
% the result is the solution NG, including the solutions of the
% polynomial.
begin scalar res,ns,testpoly,ppoly,sol,s,var,dom;
%%tr print "Intersection ";
%%tr writepri (mkquote ('list . NG),'only);
%%tr print " with ";
%%tr writepri(mkquote poly,'only);
res := ();
poly := simp poly;
var:=if not domainp numr poly then groesolmvar(numr poly,vars)
else 'constant;
loop:
if NG=() then goto finish;
ns := car NG; NG := cdr NG;
testpoly := poly;
dom := groesoldmode!* or 'RATIONAL;
groesetdmode(dom,T);
testpoly := simp prepsq testpoly;
for each u in ns do
if idp lhs u and not smemq('root_of,rhs u) then
<<s := rhs u;
testpoly := subsq(testpoly,list (lhs u . s));
>>;
groesetdmode(dom,nil);
ppoly := prepf numr testpoly;
sol := groesolvePolyv(ppoly,var);
res := append(res , for each r in sol collect append(r,ns) ) ;
goto loop;
finish:
%%tr print "Schnittresultat: ";
%%tr writeres res;
return res;
end;
symbolic procedure groesolmvar(poly,vars);
% select main variable wrt vars sequence.
<<while vars and not smember(car vars,poly) do
vars:=cdr vars;
if null vars then rerror('groebner,27,"illegal poly");
car vars>>;
% solving a single polynomial with respect to its main variable
symbolic procedure groeSolvePoly p; groesolvePolyv(p,mainvar p);
symbolic procedure groesolvePolyv(p,var);
% find the zeros for one polynomial p in the
% variable "var".
% current dmode is NIL.
(begin scalar res,u,!*convert,y,z;
if (u:=assoc(var,depl!*)) then depl!*:=delete(u,depl!*);
if !*trgroesolv then
<<writepri(" solving univariate with respect to ",'first);
writepri(mkquote var,'last);
writepri(mkquote p,'only);
>>;
for each s in groebroots!* do
if 0=reval{'difference,p,car s} then res:=cdr s;
if res then return res;
groesetdmode(groesoldmode!*,T);
u := numr simp p;
res := if !*groebnumval and UnivariatePolynomial!? u
then groeroots(p,var)
else (solveeval list (p,var))
where kord!*=nil,
alglist!* = nil . nil;
res := cdr res;
% Collect nontrivial denominator factors.
% Reorder for different local order during solveeval.
for each x in res do
<<y:=prepf (z:=reorder denr simp caddr x);
if dependsl(y,variables!*) then groepostcollectden z;
>>;
res := for each x in res collect list x;
groesetdmode(groesoldmode!*,NIL);
return res;
end) where depl!*=depl!*;
symbolic procedure UnivariatePolynomial!? fm;
domainp fm or UnivariatePolynomial!?1 (fm,mvar fm);
symbolic procedure UnivariatePolynomial!?1 (fm,v);
domainp fm or
domainp lc fm and v = mvar fm and
UnivariatePolynomial!?1(red fm,v);
symbolic procedure predecessor (r,l);
% looks for the predecessor of r in l
if not pairp l or not pairp cdr l or r = car l
then rerror(groebner,9,"No predecessor available") else
if r = cadr l then car l else predecessor (r,cdr l);
symbolic procedure ZerosetUnion(ng1,ng2);
<<%print list( "union von ",ng1,ng2;
ng1 := ZerosetUnion1(ng1,ng2);
% print list( "-->",ng1);
ng1>>;
symbolic procedure ZerosetUnion1(ng1,ng2);
% Vereinigung von Nullstellengebilden
if ng1 = () then ng2
else
if ZerosetMember(car ng1,ng2) then ZerosetUnion1(cdr ng1,ng2)
else
car ng1 . ZerosetUnion1(cdr ng1,ng2);
symbolic procedure ZerosetMember (ns,ng);
if ng = () then nil else
if ZeroEqual(ns,car ng) then ng else
ZerosetMember (ns,cdr ng);
symbolic procedure ZeroEqual(ns1,ns2);
if Zerosubset(ns1,ns2) then Zerosubset(ns2,ns1) else nil;
symbolic procedure Zerosubset(ns1,ns2);
if null ns1 then T else
if member(car ns1,ns2) then Zerosubset(cdr ns1,ns2)
else nil;
symbolic procedure groesetdmode(dmode,dir);
% Interface for switching an arbitrary domain on/off.
% Preserve complex mode. Turn on EZGCD whenever possible.
if null dmode then nil else
begin scalar !*msg,x,y;
if null dir then
<< if !*complex then y:=setdmode('complex,nil);
!*complex := nil;
if dmode!* = '!:rd!: then !*rounded :=nil;
if dmode!* then y:=setdmode(get(dmode!*,'dname),nil);
if !*groebcomplex then
<<setdmode('complex,t); !*complex:=t>>;
>>
else
<< if memq(dmode,'(complex complex!-rounded complex!-rational))
then <<!*complex := t; y:=setdmode('complex,t);
if (x:=atsoc(dmode,'((complex!-rounded . rounded)
(complex!-rational . rational)) ))
then y:=setdmode(cdr x,t)>>
else y:=setdmode(dmode,t);
if memq(dmode,'(rounded complex!-rounded)) then !*rounded :=t;
>>;
!*ezgcd := null dmode!*;
return y;
end;
symbolic procedure preduceeval pars;
% Polynomial reduction driver. u is an expression and v a list of
% expressions. Preduce calculates the polynomial u reduced wrt the list
% of expressions v.
% parameters:
% 1 expression to be reduced
% 2 polynomials or equations; base for reduction
% 3 optional: list of variables
begin scalar vars,x,u,v,w,oldorder;
scalar !*factor,!*exp,!*gsugar,!*vdpinteger;
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,3,"Empty list in Preduce");
vars := groebnervars(w,v);
if not vars then vdperr 'PREDUCE;
oldorder := vdpinit vars;
w := for each j in w collect a2vdp j;
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 := groebNormalForm(x,w, 'sort);
w := vdp2a w;
setkorder oldorder;
return if w then w else 0;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% the following code is the interface to Stan's rootfinder
%
symbolic procedure groeroots(p,x);
begin scalar r;
x := nil; r := reval{'roots,p};
% re-evaluate rhs in order to get prefix form
r := for each e in cdr r collect
list('equal,cadr e,reval caddr e);
return 'list . r;
end;
endmodule;
end;