#The following program  is written in  Maple V realease 3. 
#The following package <<wsolve>> is for solving systems 
#of polynomial equations.
#The main algorithm  is based on Wu's method. 
#Polynomial factorization over rational field is used to 
#reduce the sizes of the polynomials produced in the 
#Well-Ordering process.

#-------------------------------------------------------------
#Part I: <<wsolve>> for solving systems of polynomial equations
#-------------------------------------------------------------
readlib(factors):

# the class of poly f wrt variable ordering ord
Class := proc(f,ord)
	 local i;
	 options remember,system;
	 for i to nops(ord) do
	   if has(f,ord[i]) then RETURN(nops(ord)-i+1) fi
	 od;
	 0
	end:
# the leading variable of poly f wrt variable ordering ord
Mainvar := proc(f,ord)
	   local i;
	   options remember,system;
	   for i to nops(ord) do
	     if has(f,ord[i]) then RETURN(ord[i]) fi
	   od;
	   lprint(`Warning: lvar is called with constant`);
	   0
	  end:
# the initial of poly f wrt ord
Initial := proc(f,ord)
	   options remember,system;
	   if Class(f,ord) = 0 then 1
	   else lcoeff(f,Mainvar(f,ord)); numer("/lcoeff("))
	   fi
	   end:
# the separant of poly f wrt ord
Sept := proc(f,ord)
	options remember,system;
	if Class(f,ord) = 0 then 1
	else diff(f,Mainvar(f,ord));
	fi
	end:
# the set of all nonconstant factors of separants of polys in as
Septs := proc(as,ord)
	 local i,is,iss;
	 is := {};
	 for i in as do
	   Sept(i,ord);
	   if Class(",ord) > 0 then is := {op(is),"} fi
	 od;
	 iss := {};
	 for i in is do
	   if Class(i,ord) > 0 then iss := {op(iss),i} fi
	 od;
	 iss
	 end:
# modified rank of two polys: comparing further the rank of
# initials, the terms of initials and the terms of f and g
# when they are the same
Rank := proc(f,g,ord)
	local ind,find,cf,cg;
	options remember,system;
	find := Subrank(f,g,ord,'ind');
	if find and ind = 1 then
	  cf := nops(expand(Initial(f,ord)));
	  cg := nops(expand(Initial(g,ord)));
	  if cf < cg then true
	  elif cf = cg then
	    if nops(expand(f)) < nops(expand(g)) then 
	      true 
	    else 
	      false 
	    fi
	  else 
	    false
	  fi
	else 
	  find
	fi
	end:
# subroutine for rank
Subrank := proc(f,g,ord,ind)
	   local cf,cg;
	   options remember,system;
	   cf := Class(f,ord);
	   cg := Class(g,ord);
	   if cf = 0 then
	     if cg = 0 then ind := 1 fi; true
	   elif cf < cg then true
	   elif cf = cg then
	     cf := degree(f,Mainvar(f,ord));
	     cg := degree(g,Mainvar(g,ord));
	     if cf < cg then true
	     elif cf = cg then
	       Subrank(Initial(f,ord),Initial(g,ord),
		       ord,'ind')
	     else false
	     fi
	   else false
	   fi
	   end:
# modified pseudo division: I1^s1...Ir^sr*uu = q*vv + r,
# where I1, ..., I_r are all distinct factors of lcoeff(vv,x)
# and s1, ..., sr are chosen to be the smallest
NPrem := proc(uu,vv,x)
	 local r,v,dr,dv,l,t,lu,lv;
	 options remember,system;
	 if type(vv/x,integer) then subs(x = 0,uu)
	 else
	   r := expand(uu);
	   dr := degree(r,x);
	   v := expand(vv);
	   dv := degree(v,x);
	   if dv <= dr then 
	     l := coeff(v,x,dv); 
	     v := expand(v-l*x^dv)
	   else l := 1
	   fi;
	   while dv <= dr and r <> 0 do
	     gcd(l,coeff(r,x,dr),'lu','lv');
	     t := expand(x^(dr-dv)*v*lv);
	     if dr = 0 then r := 0 
	     else r := subs(x^dr = 0,r) 
	     fi;
	     r := expand(lu*r)-t;
	     dr := degree(r,x)
	   od;
	   r
	 fi
	 end:
# pseudo remainder of poly f wrt ascending set as
Premas := proc(f,as,ord)
	  local remd,i;
	  remd := f;
	  for i to nops(as)  do
	    remd := NPrem(remd,as[i],Mainvar(as[i],ord))
	  od;
	  if remd <> 0 then numer(remd/lcoeff(remd)) 
	  else 0 
	  fi
	  end:
##############################
#The following is for testing#
##############################
Qpremas := proc(f,as,ord)
	   local remd,i;
	   remd := f;
	   for i to nops(as) do
	     if Class(remd,ord) = Class(as[i],ord) then
	       remd := NPrem(remd,as[i],Mainvar(as[i],ord))
	     fi
	   od;
	   if remd <> 0 then numer(remd/lcoeff(remd)) 
	   else 0 
	   fi
	   end:

######################
######################
# set of non-zero remainders of polys in ps wrt ascset as
remseta := proc(ps,as,ord)
	   local i,set,r;
	   set:={};
	   for i in ps do 
	     r:=Premas(i,as,ord);
	     if r <>0  and Class(r,ord) =0
	     then RETURN({1})
	     else set:=set union {r} fi;
	   od;
	   set minus {0}
	   end:
#############################
#The following is for testing
#############################
Qremseta := proc(ps,as,ord)
	    local i,set,r;
	    set:={};
	    for i in ps do r:=Qpremas(i,as,ord);
	      if r <>0  and Class(r,ord) =0
	      then RETURN({1})
	      else set:=set union {r} fi
	    od;
	    set minus {0}
	    end:
###########################################
# search an element with lowest rank in ps
#      and assign the rest of polys to qs
###########################################
Sort := proc(ps,rank,ord,qs)
	local l,i,qs1;
	if nops(ps) = 1 then qs := []; ps[1]
	else
	  l := ps[1];
	  qs1 := [];
	  for i from 2 to nops(ps) do
	    if rank(ps[i],l,ord) then 
	      qs1 := [l,op(qs1)]; 
	      l := ps[i]
	    else qs1 := [ps[i],op(qs1)]
	    fi
	  od;
	  qs := qs1;
	  l
	fi
	end:

# the basic set of polyset ps
basset := proc(ps,ord)
	  local qs,qs1,i,b;
	  if nops(ps) < 2 then ps
	  else
	    b := Sort(ps,Rank,ord,'qs1');
	    qs := [];
	    if 0 < Class(b,ord) then
	      for i in qs1 do
		if degree(i,Mainvar(b,ord)) <
		   degree(b,Mainvar(b,ord)) 
		then
		  qs := [i,op(qs)]
		fi
	      od
	    else RETURN([b])
	    fi;
	    [op(basset(qs,ord)),b]
	  fi
	  end:

# the weak basic set of polyset ps
Wbasset := proc(ps,ord)
	   local qs,qs1,i,b;
	   if nops(ps) < 2 then ps
	   else
	     b := Sort(ps,Rank,ord,'qs1');
	     qs := [];
	     if 0 < Class(b,ord) then
	       for i in qs1 do
		 if Class(b,ord) < Class(i,ord) and
		    degree(Initial(i,ord),Mainvar(b,ord))<
		    degree(b,Mainvar(b,ord)) 
		 then
		   qs := [i,op(qs)]
		 fi
	       od
	     else RETURN([b])
	     fi;
	     [op(Wbasset(qs,ord)),b]
	   fi
	   end:

# the quasi-basic set of polyset ps
Qbasset := proc(ps,ord)
	   local qs,qs1,i,b;
	   if nops(ps) < 2 then ps
	   else
	     b := Sort(ps,Rank,ord,'qs1');
	     qs := [];
	     if 0 < Class(b,ord) then
	       for i in qs1 do
		 if Class(b,ord) < Class(i,ord) then
		   qs := [i,op(qs)]
		 fi
	       od
	     else RETURN([b])
	     fi;
	     [op(Qbasset(qs,ord)),b]
	   fi
	   end:
###########################################################
#Procedure Nrs:
#input  a polynomial set  RS:={R1,R2,...,Rn}
#output a set in which every element is a poly set. 
#         set:={set1,set2,...,sets}
###########################################################
Nrs := proc(rs,ord,test)
       local rm,rss,l1,i,j,rset;
       rss := Produce(rs,ord,test);
       if rss = {} then RETURN({}) fi;
       l1 := nops(rss);
       rset := {};
       if l1 = 1 then
	 for i in rss[1] do  rset := rset union {{i}} od; 
       RETURN(rset)
       fi;
       for i in rss[1] do  
	 for j in rss[2] do  
	   rset := rset union {{i,j}} 
	 od 
       od;
       for i from 3 to l1 do  rset := Lltl(rset,rss[i]) od;
       l1 := nops(rset);
       rm := {};
       for i to l1-1 do
	 for j from i+1 to l1 do
	   if op(i,rset) minus op(j,rset) = {} then
	      rm := rm union {op(j,rset)}
	   fi;
	   if op(j,rset) minus op(i,rset) = {} then
	      rm := rm union {op(i,rset)}
	   fi
	 od
       od;
       rset := rset minus rm;
       rset
       end:

#input  two poly set set1,set2
#output
Ltl := proc(list1,list2)
       local set,i,j;
       set := {};
       for i in list1 do  
	 for j in list2 do  
	   set := set union {{j,i}} 
	 od 
       od
       end:

Lltl := proc(llist1,list2)
      local set,i,j;
      set := {};
      for i in llist1 do
	for j in list2 do  
	  set := set union {i union {j}} 
	od
      od
      end:

Nonumber := proc(list,ord)
	    local i,ls;
	    ls := {};
	    for i in list do  
	      if 0 < Class(i,ord) then 
		 ls := ls union {i} 
	      fi 
	    od;
	    ls
	    end:

Factorlist := proc(pol)
	    local i,list1,list2;
	    list1 := factors(pol)[2];
	    list2 := {};
	    for i in list1 do  
	      list2 := list2 union {numer(i[1])} 
	    od;
	    list2
	    end:

Init := proc(p,ord)
	local pol;
	pol:=Initial(p,ord);
	if Class(pol,ord)=0 then RETURN(1) fi;
	pol
	end:

Inits := proc(bset,ord)
	 local i,list;
	 list:={};
	 for i in bset do 
	   list:=list union Factorlist(Init(i,ord)) 
	 od;
	 for i in list do 
	   if Class(i,ord)=0 then 
	      list:=list minus {i} 
	   fi 
	 od;
	 list
	 end:

Inits1 := proc(bset,ord)
	  local i,list;
	  list:={};
	  for i  to nops(bset)-1 do 
	    if Class(Init(i,ord),ord)>1 then 
	       list:=list union Factorlist(Init(i,ord)) 
	    fi 
	  od;
	  for i in list do 
	    if Class(i,ord)=0 then 
	       list:=list minus {i} 
	    fi 
	  od;
	  list
	  end:

nums := proc(ps,ord)
	local i,n;
	n:=0;
	for i in ps do if Class(i,ord)=1 then n:=n+1 fi od;
	n
	end:

Produce := proc(rs,ord,test)
	   global remfac;
	   local i,j,list1,list2;
	   list2 := {};
	   for i in rs do
	     list1 := Nonumber(Factorlist(i),ord);
	     remfac := remfac union (list1 intersect test);
	     list1 := list1 minus test;
	     for j in list1 do
	       if Class(j,ord) = 0 then 
		  list1 := list1 minus {j} 
	       fi
	     od;
	     list2 := list2 union {list1}
	   od;
	   for j in list2 do  
	     if j = {} then RETURN({}) fi 
	   od;
	   list2
	   end:

Reverse := proc(list)
	   local i,n,list1;
	   n := nops(list); 
	   list1 := [seq(list[n-i+1],i = 1 .. n)]; 
	   list1
	   end:

NTest := proc(ps,as,ord)
	 local i;
	 for i in ps do  
	   if normalf(i,as,ord) = 0 then 
	      RETURN(0) 
	   fi 
	 od;
	 1
	 end:

TestQ := proc(ps,as,ord)
	 global remfac;
	 local i;
	 for i in ps do  
	   if normalf(i,as,ord) = 0 then 
	      remfac:=remfac union {i};
	      RETURN(0) 
	   fi 
	 od;
	 1
	 end:

Cs1 := proc(ps,ord,nzero,test,T)
       global AS,iindex;
       local asc,i,j,rm,cset,test1,ps1,bset,rset;
       options remember,system;
       if nums(ps,ord)>1 then RETURN({}) fi;
       rm := {};
       cset := {};
       if {op(ps)} intersect nzero <> {} then RETURN({}) fi;
       if nops(nzero) <> 0 then 
	  if NTest(nzero,ps,ord) = 0 then 
	     print(0);
	     RETURN({}) 
	  fi  
       fi;
       if nops(test)  <> 0 then 
	  if TestQ(test ,ps,ord) = 0 then 
	     print(0);
	     RETURN({}) 
	  fi  
       fi;
       bset := T.`basset`(ps,ord);
       rset := T.`remseta`([op(ps)],bset,ord);
       if rset={1} then Return({0}) fi;
       if rset = {} then
	 iindex := iindex+1;
	 AS[iindex] := bset;
	 print(`A New Component`);
	 RETURN({bset})
       fi;
       test1 := test union Inits(bset,ord);
       rset := Nrs(rset,ord,nzero union test1);
       if rset = {} then RETURN({}) fi;
       cset:=map(Cs1,map(`union`,rset,{op(bset)}),
		 ord,nzero,test1,T);
       map(op,cset);
       end:

CS1 := proc(ps,ord,nzero,T)
       local  pset,cset,order,nonzero,asc;
       options remember, system;
       if nargs < 1 then ERROR(`too few arguments`)
       elif nargs>4 then ERROR(`too many arguments`)
       fi;
       if nargs<2 then order:=[op(indets(ps))] 
       else order:=[op(ord)] 
       fi;
       if nargs<3 then nonzero:={} else nonzero:=nzero fi;
       if nargs<4 then asc:=`` else asc:=T fi;
       cset:={};
       pset:=Nrs(ps,order,nonzero);
       cset:=map(Cs1,pset,order,nonzero,{},asc);
       [op(map(op,cset))];
       end:

Cs2 := proc(ps,ord,nzero,test,T)
       local cset,cset1,i,j,rs,rs1;
       options remember,system;
       if nums(ps,ord)>1 then RETURN({}) fi;
       rs1:={};
       cset := Cs1(ps,ord,nzero,test,T);
       cset1:=cset;
       for i in cset1 do rs:=remseta([op(ps)],i,ord);
	   if rs<>{} then 
	      rs1:=rs1 union {rs union {op(i)}};
	      cset:=cset minus {i} 
	   fi
       od;
       if rs1={} then RETURN(cset) fi;
       for j in rs1 do  
	 cset:=cset union Cs2(ps union j,ord,nzero,test,T) 
       od;
       cset
       end:

Cs := proc(ps,ord,nzero,T)
      local cset,remf,ps1,i,nzeros;
      global remfac;
      options remember,system;
      if nums(ps,ord)>1 then RETURN({}) fi;
      remfac:={};
      cset:=Cs2({op(ps)},ord,{op(nzero)},{},T);
      remf:=remfac minus {op(nzero)};
      for i to nops(remf) do
	ps1 := {op(ps)} union {remf[i]};
	if i = 1 then nzeros := {op(nzero)} 
	else nzeros := nzeros union {remf[i-1]} 
	fi;
	cset := cset union Cs(ps1,ord,nzeros,T)
      od;
      cset
      end:

CS := proc(ps,ord,nzero,T)
      local i,pset,cset,order,nonzero,asc;
      options remember, system;
      if nargs < 1 then ERROR(`too few arguments`)
      elif nargs > 4 then ERROR(`too many arguments`)
      fi;
      if nargs < 2 then order:=[op(indets(ps))] 
      else order:=[op(ord)] 
      fi;
      if nargs < 3 then nonzero:={} else nonzero:=nzero fi;
      if nargs < 4 then asc:=`` else asc:=T fi;
      cset:={};
      pset:=Nrs(ps,order,nonzero);
      for i to nops(pset) do 
	  cset:=cset union Cs(pset[i],order,nonzero,asc) 
      od;
      [op(cset)]
      end:

Css := proc(ps,ord,nzero,T)
       global remfac;
       local cset1,cset,i,j,nzero1, ps1,ps2,Is,Ninits;
       options remember,system;
       remfac := {};
       cset1 := Cs(ps,ord,nzero,T);
       cset := cset1;
       for i in cset1 do
	 if Class(i[1],ord)=1 then 
	    Is:=Inits1(i,ord) 
	 else Is := Inits(i,ord) 
	 fi;
	 Ninits := Is minus {op(nzero)};
	 for j to nops(Ninits) do
	   ps1 := ({op(ps)} union {op(i)}) union {Ninits[j]};
	   if j = 1 then nzero1 := {op(nzero)}
	   else nzero1 := nzero1 union {Ninits[j-1]}
	   fi;
	   cset := cset union Css(ps1,ord,nzero1,T)
	 od
       od;
       cset
       end:

CSS := proc(ps,ord,nzero,T)
       local asc,i,pset,cset,nonzero,order;
       options remember,system;
       if nargs < 2 then order:=[op(indets(ps))] 
       else order:=[op(ord)] 
       fi;
       if nargs < 3 then nonzero:={} 
       else nonzero:=nzero 
       fi;
       if nargs < 4 then asc:=`` else asc:=T fi;
       cset:={};
       pset:=Nrs(ps,order,nonzero);
       for i to nops(pset) do 
	   cset:=cset union Css(pset[i],order,nonzero,asc) 
       od;
       cset
       end:

Css1 := proc(ps,ord,test)
	local  cset,i,j,septs,cset1;
	cset:=CSS(ps,ord,test);
	septs:={};
	cset1:=cset;
	for i in cset1 do 
	    septs:=Septs(i,ord) minus Inits(i,ord);
	    if septs<>{} then
	       for j in septs do 
		   cset:=cset union Css1(ps union 
			 {op(i)} union {j},ord,test) 
	       od
	     fi
	od;
	cset
	end:

#Gsolve:=proc(ps,ord,test)
#	gsolve(ps,ord,test)
#	end:

wsolve:=proc() 
	[op(CSS(args))] 
	end:

e_val:=proc(ps,ord,test)
       Css1(ps,ord,test)
       end:

iindex:=0:
remfac:={}:

#---------------------------------------------------
#Part II: solving the prepared the system
#---------------------------------------------------

#solve1 for solving a system of tri set about the mainvar;
solve1 := proc(cs,ord)
	  local sol,i;
	  sol:=[];
	  for i from  nops(cs) by -1  to 1 do
	    sol:=[op(sol),[solve(cs[i],
		  {Mainvar(cs[i],ord)})]]
	  od;
	  ltime(sol)
	  end:

Ltl_ := proc(list1,list2)
	local set,i,j;
	set := {};
	for i in list1 do  
	    for j in list2 do  
		set := set union {[op(i),op(j)]} 
	    od 
	od
	end:

Lltl_ := proc(llist1,list2)
	 local set,i,j;
	 set := {};
	 for i in llist1 do
	     for j in list2 do  
	       set := set union {[op(i),op(j)]} 
	     od
	  od
	  end:
 
sol_ := proc(sol)
	local list1,i;
	list1:=sol;
	for i from nops(list1) by -1 to 2 do
	  list1:=subs(list1[i]=subs(op(Reverse(
		 [list1[1..i-1]])),list1[i] ),list1)
	od;
	simplify(list1);
	end:

ltime := proc(pss)
	 local i,l,ls;
	 l:=nops(pss);
	 ls:=Ltl_(pss[1],pss[2]);
	 if l > 2 then
	    for i from 3 to l do
	      ls:=Lltl_(ls,pss[i])
	    od;
	 fi;
	 ls
	 end:

w_eval := proc(ps,sol)
	  subs(op(sol),ps);
	  end:

#simplifyas := proc(as,ord,nzero)
#	      local cs;
#	      if type(as[1],list) then 
#		 cs:=map(simplifyas,as,ord,nzero)
#	      else
#		 print(`Simplify the ascending set`);
#		 cs:=gsolve(as,ord,nzero)
#	      fi;
#	      cs
#	      end:

solveas := proc(as,ord)
	   local soln;
	   if type(as[1],list) then 
	      soln:=map(solveas,as,ord)
	   else
	      soln:=map(sol_,solve1(as,ord))
	   fi;
	   soln;
	   end:

wsolve_s := proc(ps,ord,nzero)
	    local css,soln;
	    print(`Step 1:Compute the ascending sets`);
	    css:=wsolve(ps,ord,nzero);
	    print(`Step 2:Compute the solutions`);
	    soln:=solveas(css,ord);
	    soln
	    end:

#-----------wsolve help ------

`help/text/wsolve` := TEXT(
`FUNCTION: wsolve - prepare the given algebraic system for solving`,
`   `,
`CALLING SEQUENCE:`,
`   wsolve(F)`,
`   wsolve(F, X)`,
`   wsolve(F, X,nonzero)`,
`   `,
`PARAMETERS:`,
`   F       - set or list of polynomials in indeterminates X`,
`   X       - list of indeterminates (not including parameters)`,
`   nonzero - set of polynomials in indeterminates X`,
`   `,
`SYNOPSIS:`,
`- The command wsolve(F,X,nonzero) computes a collection of `,
`ascending sets  corresponding to F.`,
`   `,
`- First, the system corresponding to F is subdivided by `,
`factorization.`,
`   `,
`- Then, each subsystem is passed to a variant of Wu's `,
`algorithm which factors all intermediate results.`,
`  `,
`- The result is a set of subsystems whose roots which are `,
`not the zeros of`,
`the initials (leading coefficients to the main variable) `,
`are those of the original system, but whose variables have`, 
`been successively eliminated`,
`  `,
`- If desired, the solveas(as,ord) function may then be applied`, 
`to solve the subsystems`,
`  `,
`- If desired, the simplifyas(as,ord,nonzero) function may `,
`then be applied to simplify the ascending set.`,
`  `,
`- The set nonzero may be used to prevent certain quantities `,
`from being considered in roots; however, this may not stop `,
`all such solutions.`,
`  `,
`  `,
`EXAMPLES:`,
`> F := [x^2 - 2*x*z + 5, x*y^2 + y*z^3, 3*y^2 - 8*z^3]:`,
`> wsolve(F,[z,y,x]);`,
`                                                 `,
`              2                `,
` {[z, y, x  + 5],          `,
`                               `,
`   2                            6       5       4        2       `,
` [x  - 2 x z + 5, 3 y + 8 x, 3 x  - 64 x  + 45 x  + 225 x  + 375]}`,
`       `
):