% Supplement to the paper
% About the Polynomial System Solve Facility of Axiom, Macsyma,
% Maple, Mathematica, MuPAD, and Reduce
% H.-G. Graebe | Univ. Leipzig


% ####################################
%	Some useful code
% ####################################

Axiom:

)set message time on
numsolve1(polys) == complexSolve(polys,0.0001); 
numsolve2(sol) == concat([complexSolve(u,0.0001) for u in sol]);
resubst(sol,polys) == [[subst(u,v) for u in polys] for v in sol];

Macsyma:

showtime:true$
display2d:false$

resubst(sol,polys):=
map(lambda([elem],ratsimp(subst(elem,polys))),sol)$


Maple:

tt:=time();
sol:=[solve(polys,vars)];
time()-tt;
nops(sol);
lprint(sol);

expandsol := proc(sol) map(allvalues,sol,dependent) end;
numsolve2 := proc(sol) map(evalf,expandsol(sol)) end;
resubst:=proc(sol,polys) map(simplify@subs,sol,polys) end;

Mathematica:

eqn[n_List]:=(n==Table[0,{Length[n]}])
time:=res[[1]]
numsolve1[polys_]:=NSolve[eqn[polys],vars]//Timing
numsolve2[sol_]:=sol//N//Timing
resubst[sol_,polys_]:=polys/.sol//Simplify
res=Solve[eqn[polys],vars]//Timing
sol=res[[2]]
sol//InputForm

procedure sol2eqns(sol);
  for each a in sol collect for each b in a collect (lhs b - rhs b);

MuPAD:

tt:=time((sol:=solve(polys,vars))); sol;
tt:=time((sol1:=numsolve1(polys,vars))); nops(sol1);
tt:=time((sol2:=map(sol,op@allvalues))); nops(sol2);

numsolve1:=proc(polys,vars) begin float(hold(solve)(polys,vars))
end_proc;

mysubs:=proc(a,b) begin simplify(subs(b,op(a))) end_proc;
resubst:=proc(sol,polys) begin map(sol,mysubs,polys) end_proc;

Reduce:

off nat; on time;

procedure sol2eqns(sol); 
  for each a in sol collect for each b in a collect (lhs b - rhs b);

procedure numsolve1(); 
	<< on rounded; write solve(polys,vars); off rounded; >>;

procedure numsolve2(); 
	<< on rounded; 
	write for each x in expand_cases(sol) collect sub(x,x);
	off rounded; >>;

procedure resubst(sol,polys); 
	for each u in sol collect sub(u,sub(u,polys));


% ####################################
%	A first example (Arn 4)
% ####################################

vars:={w,x,y,z};

polys:={w + x + y + z,
w*z + w*x + x*y + y*z,
x*w*z + w*y*z + w*x*y +
x*y*z,w*x*y*z - 1};

% ####################################
%	zero dimensional examples
% ####################################

% Ex 1:

vars:={x,y,z};

polys:={x^2 + y + z - 3,x + y^2 + z - 3,x + y + z^2 - 3};

% -----------

% Ex 2:

vars:={x,y};

polys:={x^4+x+1,y^4+x+1};

% -----------

% Ex 3:

vars:={u0,u1,u2,u3};

polys:={u0+2*u1+2*u2+2*u3-1, u0^2+2*u1^2+2*u2^2+2*u3^2-u0,
2*u0*u1+2*u1*u2+2*u2*u3-u1, 2*u0*u2+u1^2+2*u1*u3-u2};

% -----------

% Ex 4:

vars:={x,y,z};

polys:={x^3 + y + z - 3,x + y^3 + z - 3,x + y + z^3 - 3};

% -----------

% Ex 5:

vars:={v, w, x, y, z};

polys:={v + w + x + y + z, 
v*w + v*z + w*x + x*y + y*z,
v*w*x + v*w*z + v*y*z + w*x*y + x*y*z,
v*w*x*y + v*w*x*z + v*w*y*z + v*x*y*z + w*x*y*z,
v*w*x*y*z - 1};

% -----------

% Ex 6:

vars:={x,y,z};

polys:={x^3 + y^2 + z - 3,x + y^3 + z^2 - 3,x^2 + y + z^3 - 3};

% -----------

% Ex 7:

vars:={u0,u1,u2,u3,u4};

polys:={u0+2*u1+2*u2+2*u3+2*u4-1,
u0^2+2*u1^2+2*u2^2+2*u3^2+2*u4^2-u0,
2*u0*u1+2*u1*u2+2*u2*u3+2*u3*u4-u1,
2*u0*u2+u1^2+2*u1*u3+2*u2*u4-u2,
2*u0*u3+2*u1*u2+2*u1*u4-u3};

% #################################################
%	zero dimensional examples with parameters
% #################################################

vars:={x,y};

polys:={a*x+y-1, x+a*y-1};

% ----------- parametric A4, pars:={z};

% Ex 8:

vars:={w,x,y}; 

polys:={w + x + y + z,
w*z + w*x + x*y + y*z,
x*w*z + w*y*z + w*x*y + x*y*z,
w*x*y*z - 1};

% ----------- [BGK,example 6], pars:={a,u,v,w};

% Ex 9:

vars:={x,y,z,t}; 

polys:= {t-(a-v), x+y+z+t-(u+w+a), x*z+x*t+y*z+z*t-(u*a+u*w+w*a),
x*z*t-(u*w*a)};

% ----------- The ROMIN 3R-robot [glopez,recio], pars:={a,b,c,l1,l2};

% Ex 10:

vars:={s1,s2,s3,c1,c2,c3,d}; 

polys:={a + d*s1, - b + c1*d, c2*l2 + c3*l3 - d, - c + l2*s2 + l3*s3,
c1^2 + s1^2 - 1, c2^2 + s2^2 - 1, c3^2 + s3^2 - 1};


% ##############################
%	Arbitrary systems
% ##############################

% Ex E1: Computing the extrema of f(x,y)=x^3*y^2*(6-x-y)

vars:={x,y};

m := {x^2*y^2*( - 4*x - 3*y + 18), x^3*y*( - 2*x - 3*y + 12)};

% ----------- 

% Ex G1: Gerdt et al. [6, (4)] : Seventh order KdV type equation.

vars:= {l7, l6, l5, l4, l3, l2, l1};

polys:= {l1*(2*l4 - l5 + 2*l6), 
 - 20*l1^3 + 10*l1^2*l2 - 2*l1^2*l3 + 70*l1*l4 - 35*l2*l4 + 7*l3*l4,
 6*l1^2*l4 - 2*l1^2*l5 + 2*l1^2*l6 - 21*l4^2 + 7*l4*l5 - 7*l4*l6,
 15*l1^3 - 7*l1^2*l2 - 6*l1^2*l3 + 5*l1*l2^2 + 4*l1*l2*l3 - 42*l1*l4 +
21* l1*l5 - 63*l1*l6 - 2*l2^3 - 14*l2*l5 + 42*l2*l6 + 147*l7, 
 - 135*l1^4 + 45*l1^3*l2 - 9*l1^3*l3 + 602*l1^2*l4 + 28*l1^2*l5 -
196*l1*l2*l4 - 14*l1*l2*l5 + 70*l1*l3*l4 - 28*l1*l3*l5 - 2205*l1*l7 -
14*l2^2*l4 + 14*l2^2*l5 + 735*l2*l7 - 147*l3*l7 - 98*l4*l5 + 294*l4*l6
+ 98*l5^2 - 294*l5*l6,
 36*l1^3*l4 - 9*l1^3*l5 + 6*l1^3*l6 - 28*l1^2*l7 + 14*l1*l2*l7 +
28*l1*l3*l7 - 168*l1*l4^2 + 42*l1*l4*l5 - 28*l1*l4*l6 - 14*l2^2*l7 +
588*l4*l7 - 245*l5*l7 + 392*l6*l7,
 l1*( - 5*l1^2 + 13*l1*l2 - l1*l3 - 6*l2^2 + 2*l2*l3 - 28*l4 + 14*l6),
 l1*(5*l1*l3 - 3*l2*l3 + l3^2 - 28*l4 + 14*l6),
 2*l1*(11*l1*l4 - 5*l1*l5 - 8*l1*l6 - l2*l4 + 3*l2*l5 + 2*l2*l6 -
l3*l4 - l3*l5 + 21*l7),
 4*l1*( - 33*l1*l4 + 10*l1*l5 + 39*l1*l6 + 17*l2*l4 - 6*l2*l5 -
22*l2*l6 - 5*l3*l4 + 2*l3*l5 + 7*l3*l6 - 63*l7),
 l1*(15*l1*l7 - 30*l2*l7 + 12*l3*l7 - 20*l4^2 + 8*l4*l5 + 10*l4*l6 -
4*l5*l6),
 l1*( - 6*l1*l7 + 12*l2*l7 - 6*l3*l7 + 4*l4^2 - 2*l4*l5 + 2*l4*l6 +
l5*l6 - 2*l6^2),
 4*l1*l7*(2*l4 - l5 + 2*l6)};

% ----------- 

% Ex G6: Gerdt et al. [6, (8)] : Two coupled non linear eq's of KdV type.
%       Dort t=L1 ,x=L5 ,y=L3 ,z=L4 

vars:={t,x,y,z};

polys:={
 -2*z^3*t+(3*z^2*t-2*z^2-6*z*y*t+6*z*y+6*y^2*t-6*y^2)*x-z*t*x^2,
18*z^3*t^2-9*z^3*t-18*z^2*y*t^2+18*z^2*y*t+18*z*y^2*t^2-18*z*y^2*t+
        (-27*z^2*t^2+24*z^2*t-5*z^2+63*z*y*t^2-78*z*y*t+15*z*y-63*y^2*t^2+
        78*y^2*t-15*y^2)*x+9*z*t^2*x^2,
-8*z^4*t+(6*z^3*t-6*z^3-12*z^2*y*t+12*z^2*y+12*z*y^2*t-12*z*y^2)*x+
        (5*z^2*t-4*z^2-18*z*y*t+18*z*y+18*y^2*t-18*y^2)*x^2-3*z*t*x^3,
(3t-5)*z^2*y-15*(t-1)*z*y^2+10*(t-1)*y^3+(z*y+3*y^2*t-3*y^2)*x-y*t*x^2};

% ----------- 

% Ex Go: Gonnet, ACM SIGSAM Bulletin 17 (1983), 48 - 49
%	homogenized with a1 ( original : a1=1)

vars:={a0,a1,a2,a3,a4,a5,b0,b1,b2,b3,b4,b5,c0,c1,c2,c3,c4,c5};

polys:={a4*b4, a5*b1+a1*b5+a4*b3+a3*b4, a2*b2,a5*b5,
(a0+a1+a4)*b2+a2*(b0+b1+b4)+c2,
(a0+a1+a4)*(b0+b1+b4)+(a3+a5)*b2+a2*(b3+b5)+c0+c1+c4,
(a3+a5)*(b0+b1+b4)+(b3+b5)*(a0+a1+a4)+c3+c5-a1^2,
(a3+a5)*(b3+b5), a5*(b3+b5)+b5*(a3+a5),
b5*(a0+a1+2*a4)+a5*(b0+b1+2*b4)+a3*b4+a4*b3+c5,
a4*(b0+b1+2*b4)+a2*b5+a5*b2+(a0+a1)*b4+c4,
a2*b4+a4*b2, a4*b5+a5*b4, 2*a3*b3+a3*b5+a5*b3,
c3+b3*(a0+2*a1+a4)+a3*(b0+2*b1+b4)+a1*b5+a5*b1,
c1+(a0+2*a1+a4)*b1+a2*b3+a3*b2+a1*(b0+b4),
a2*b1+a1*b2, a5*b3+a3*b5, a1*b4+a4*b1};

% ----------- 

% Ex Kri:

comment 

Subject: Polynomial Equations 
Date: Sat, 01 Jun 1996 23:09:19 -0700 
From: "Elizabeth S. Krider" <ekrider@cco.caltech.edu>
	241, S., Wilson, Ave, #103, Pasadena, CA, 91106
Organization: California Institute of Technology, Pasadena 
Newsgroups: sci.math.symbolic

end comment;

PM0 := -6*c*g^2*u*p*v^5-2*(2*c*d+a+b*d+b+c*d^2+w*b*f-w+2
*w*c*d*f+2*w*c*f+w^2*b*g+2*w^2*c*g+2*w^2*c*d*g+w^2*c*f^2+2*w
^3*c*f*g+c*g^2*w^4)*u*p*v-2*(b*g+2*c*g+2*c*d*g+c*f^2+6*w
*c*f*g+6*c*g^2*w^2)*u*p*v^3-2*c*u^3*p*v;

PM1 := -6*c*g*(2*f+5*g*w)*u*p*v^5-2*(2*c*d+a+b*d+b+c*d^2+
w*b*f-w+2*w*c*d*f+2*w*c*f+w^2*b*g+2*w^2*c*g+2*w^2*c*d*g+w^2
*c*f^2+2*w^3*c*f*g+c*g^2*w^4)*w*u*p*v-2*(b*f-1+2*c*d*f+
2*c*f+3*w*b*g+6*w*c*g+6*w*c*d*g+3*w*c*f^2+12*w^2*c*f*g+10*c*g^
2*w^3)*u*p*v^3-2*w*c*u^3*p*v;

PM2 := -6*(b*g+2*c*g+2*c*d*g+c*f^2+10*w*c*f*g+15*c*g^2*w^2)*
u*p*v^5-2*(2*c*d+a+b*d+b+c*d^2+w*b*f-w+2*w*c*d*f+2*w*c*f+w
^2*b*g+2*w^2*c*g+2*w^2*c*d*g+w^2*c*f^2+2*w^3*c*f*g+c*g^2*w^4)*
w^2*u*p*v-2*(a+b-3*w+6*w*c*f+15*c*g^2*w^4+6*w^2*c*f^2
+12*w^2*c*g+6*w^2*b*g+3*w*b*f+12*w^2*c*d*g+20*w^3*c*f*g+6*w*c
*d*f+b*d+c*d^2+2*c*d)*u*p*v^3-30*c*g^2*u*p*v^7-2*
w^2*c*u^3*p*v-2*c*u^3*p*v^3;

ME0 := -6*c*g^2*u*p*v^5-2*(a+b*d+c*d^2-w+w*b*f+2*w*c*d*
f+w^2*b*g+2*w^2*c*d*g+w^2*c*f^2+2*w^3*c*f*g+c*g^2*w^4)*u*
p*v-2*(b*g+2*c*d*g+c*f^2+6*w*c*f*g+6*c*g^2*w^2)*u*p*
v^3-2*c*u^3*p*v;

ME1 := -6*g*(b*g+3*c*d*g+3*c*f^2+15*w*c*f*g+15*c*g^2*w^2)*u*
p*v^5-2*(a+b*d+c*d^2-w+w*b*f+2*w*c*d*f+w^2*b*g+2*w^2*c*d*
g+w^2*c*f^2+2*w^3*c*f*g+c*g^2*w^4)*(d+f*w+g*w^2)*u*p*
v-2*(-f-3*g*w+3*c*d*f^2+3*w*c*f^3+18*w^2*c*d*g^2+30*w^3*c*f
*g^2+6*g*w*b*f+18*g*w*c*d*f+18*g*w^2*c*f^2+g*a+b*f^2+6*w^2*b*g^
2+15*c*g^3*w^4+2*g*b*d+3*g*c*d^2)*u*p*v^3-30*c*g^3*u*
p*v^7-2*(b+3*c*d+3*w*c*f+3*w^2*c*g)*u^3*p*v-6*c*g*
u^3*p*v^3;

ME2 := -2*(a+b*d+c*d^2-w+w*b*f+2*w*c*d*f+w^2*b*g+2*w^2*c*d*g+
w^2*c*f^2+2*w^3*c*f*g+c*g^2*w^4)*(d+f*w+g*w^2)^2*u*p*
v-2*(-6*w*g*d+120*w^3*c*d*f*g^2+12*w*c*d*f^3+60*w^4*c*d*g^3
+40*w^3*c*f^3*g+90*w^4*c*f^2*g^2+84*w^5*c*f*g^3+6*a*g*f*w+18*b*
d*g*f*w+18*b*d*g^2*w^2+36*c*d^2*g*f*w+36*c*d^2*g^2*w^2+18*w^2*
b*f^2*g+30*w^3*b*f*g^2+72*w^2*c*d*f^2*g+a*f^2-10*g^2*w^3-3*w*f^
2-2*d*f+3*b*d^2*g+3*b*d*f^2+4*c*d^3*g+6*c*d^2*f^2-12*g*f*w^2+3*w*
b*f^3+15*w^4*b*g^3+6*w^2*c*f^4+28*c*g^4*w^6+2*a*g*d+6*a*g^2*w^2
)*u*p*v^3-6*(12*g*c*d*f^2+20*g*w*c*f^3+15*g^2*w*b*f+60*g
^2*w*c*d*f+60*g^3*w^2*c*d+90*g^2*w^2*c*f^2+140*g^3*w^3*c*f+g^2*
a-5*g^2*w-2*g*f+c*f^4+3*g^2*b*d+6*g^2*c*d^2+15*g^3*w^2*b+70*g^4*c
*w^4+3*g*b*f^2)*u*p*v^5-210*c*g^4*u*p*v^9-30*g^2*
(b*g+4*c*d*g+6*c*f^2+28*w*c*f*g+28*c*g^2*w^2)*u*p*v^7-
2*(a+3*b*d+6*c*d^2-w+3*w*b*f+12*w*c*d*f+3*w^2*b*g+12*w^2*c*d*g+
6*w^2*c*f^2+12*w^3*c*f*g+6*c*g^2*w^4)*u^3*p*v-6*(b*g+4*
c*d*g+2*c*f^2+12*w*c*f*g+12*c*g^2*w^2)*u^3*p*v^3-36*c*g^
2*u^3*p*v^5-6*c*u^5*p*v;

vars:={a,b,c,d,f,g,u,v,w,p};

polys:={PM0,PM1,PM2,ME0,ME1,ME2};

% ----------- 

% Ex G7: Gerdt, Zharkov [2]. The only ex ( eq's (6) )

vars:={a2,b2,a4,b4,a1,b1,a3,b3,a0,b0};

polys:={
 a1*a3 - a1*a4 - a4*b3 + a4*b4,
 2*a0*a1^2*a3 - a0*a1^2*a4 - a0*a1*a2*b2 - 6*a0*a2*b2*b3 +
12*a0*a3^2*b2 - 4*a0*a3*a4*b2 - a0*a4^2*b2 - a0*a4*b1*b2 -
2*a1^2*a3*b0 + a1^2*a4*b0 + a1*a2*b0*b2 - 2*a3*a4*b0*b2 + a4^2*b0*b2
+ a4*b0*b1*b2,
 a0*a1^2*a2 - 2*a0*a1*a2*b3 + a0*a1*a2*b4 + 6*a0*a2*a3*b2 +
a0*a2*a4*b2 - 12*a0*a2*b3^2 + 6*a0*a2*b3*b4 - 2*a0*a4*b1*b3 +
a0*a4*b1*b4 - a1^2*a2*b0 + 2*a1*a2*b0*b3 - a1*a2*b0*b4 - a2*a4*b0*b2
+ 2*a4*b0*b1*b3 - a4*b0*b1*b4,
 a0*a1*a4 + 3*a0*a2*b2 + 3*a0*a3*b3 + a0*a4*b3 - a1*a4*b0 - a4*b0*b3,
 5*a0^2*a1*a3 - 5*a0^2*a1*a4 - 6*a0^2*a2*b2 + 4*a0^2*a3*b3 +
8*a0^2*a3*b4 - 7*a0^2*a4*b3 + a0^2*a4*b4 - 10*a0*a1*a3*b0 +
10*a0*a1*a4*b0 - 12*a0*a2*b0*b2 + 16*a0*a3*b0*b3 - 10*a0*a3*b0*b4 +
8*a0*a4*b0*b3 - 2*a0*a4*b0*b4 + 5*a1*a3*b0^2 - 5*a1*a4*b0^2 -
2*a3*b0^2*b3 + 2*a3*b0^2*b4 - a4*b0^2*b3 + a4*b0^2*b4,
 3*a0^4*a1*a3 - 4*a0^4*a1*a4 - 6*a0^4*a2*b2 - 6*a0^4*a3*b3 -
9*a0^4*a3* b4 - 2*a0^4*a4*b3 + 9*a0^4*a4*b4 - 9*a0^3*a1*a3*b0 +
10*a0^3*a1*a4*b0 - 45*a0^3*a2*b0*b2 - 45*a0^3*a3*b0*b3 -
9*a0^3*a3*b0*b4 + 32*a0^3*a4*b0*b3 - 18*a0^3*a4*b0*b4 +
9*a0^2*a1*a3*b0^2 - 6*a0^2*a1*a4*b0^2 - 27*a0^2*a2*b0^2*b2 -
27*a0^2*a3*b0^2*b3 + 18*a0^2*a3*b0^2*b4 - 30*a0^2*a4*b0^2*b3 +
9*a0^2*a4*b0^2*b4 - 3*a0*a1*a3*b0^3 - 2*a0*a1*a4*b0^3 -
3*a0*a2*b0^3*b2 - 3*a0*a3*b0^3*b3 - a0*a4*b0^3*b3 + 2*a1*a4*b0^4 +
a4*b0^4*b3,
 - a3*b4 + a4*b4 + b1*b3 - b1*b4,
 a0*a1*a2*b4 + a0*a2*b1*b2 - 2*a0*a2*b3*b4 + a0*a2*b4^2 - 2*a0*b1^2*b3
+ a0*b1^2*b4 - a1*a2*b0*b4 - 6*a2*a3*b0*b2 - a2*b0*b1*b2 +
12*a2*b0*b3^2 - 4*a2*b0*b3*b4 - a2*b0*b4^2 + 2*b0*b1^2*b3 -
b0*b1^2*b4,
 2*a0*a1*a3*b4 - a0*a1*a4*b4 - a0*a2*b2*b4 + 2*a0*a3*b1*b2 -
a0*a4*b1*b2 - a0*b1^2*b2 - 2*a1*a3*b0*b4 + a1*a4*b0*b4 + 6*a2*b0*b2*b3
+ a2*b0*b2*b4 - 12*a3^2*b0*b2 + 6*a3*a4*b0*b2 - 2*a3*b0*b1*b2 +
a4*b0*b1*b2 + b0*b1^2*b2,
 - a0*a3*b4 - a0*b1*b4 + 3*a2*b0*b2 + 3*a3*b0*b3 + a3*b0*b4 +
b0*b1*b4,
 - 2*a0^2*a3*b3 - a0^2*a3*b4 + 2*a0^2*a4*b3 + a0^2*a4*b4 +
5*a0^2*b1*b3 - 5*a0^2*b1*b4 - 12*a0*a2*b0*b2 + 16*a0*a3*b0*b3 +
8*a0*a3*b0*b4 - 10*a0*a4*b0*b3 - 2*a0*a4*b0*b4 - 10*a0*b0*b1*b3 +
10*a0*b0*b1*b4 - 6*a2*b0^2*b2 + 4*a3*b0^2*b3 - 7*a3*b0^2*b4 +
8*a4*b0^2*b3 + a4*b0^2*b4 + 5*b0^2*b1*b3 - 5*b0^2*b1*b4,
 a0^4*a3*b4 + 2*a0^4*b1*b4 - 3*a0^3*a2*b0*b2 - 3*a0^3*a3*b0*b3 -
a0^3*a3*b0*b4 - 3*a0^3*b0*b1*b3 - 2*a0^3*b0*b1*b4 - 27*a0^2*a2*b0^2*b2
- 27*a0^2*a3*b0^2*b3 - 30*a0^2*a3*b0^2*b4 + 18*a0^2*a4*b0^2*b3 +
9*a0^2*a4*b0^2*b4 + 9*a0^2*b0^2*b1*b3 - 6*a0^2*b0^2*b1*b4 -
45*a0*a2*b0^3*b2 - 45*a0*a3*b0^3*b3 + 32*a0*a3*b0^3*b4 -
9*a0*a4*b0^3*b3 - 18*a0*a4*b0^3*b4 - 9*a0*b0^3*b1*b3 +
10*a0*b0^3*b1*b4 - 6*a2*b0^4*b2 - 6*a3*b0^4*b3 - 2*a3*b0^4*b4 -
9*a4*b0^4*b3 + 9*a4*b0^4*b4 + 3*b0^4*b1*b3 - 4*b0^4*b1*b4};