#This version is modified by Zhijie Chen in May 2003 for Maple v4
#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
# 输出 f 的类
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
# 输出 f 的主变量
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
# f 的初式
Initial := proc(f,ord)
	   options remember,system;
	   if Class(f,ord) = 0 then 1
	   else lcoeff(f,Mainvar(f,ord)); numer(%/lcoeff(%))
	   fi
	   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
# 比较 f,g 的序(类小的居前), f<g 时返回 true, 否则 false
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
# 同上比较 f,g, f<g 返回 true, 但无法确定时 ind=1
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
# 多项式 f 关于升列 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:

# set of non-zero remainders of polys in ps wrt ascset as
# 多项式组 ps 关于升列 as 的余式组, 输出 {1} 表示矛盾方程组
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:
###########################################
# search an element with lowest rank in ps
#      and assign the rest of polys to qs
###########################################
# 输出多项式组 ps 中序最小的多项式, 剩下的多项式放在 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
# 生成多项式组 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:


###########################################################
#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}
###########################################################
# 把多项式组的各多项式因式分解, 除去 test 中的因子, 然后分解成一系列
# 多项式组, 使得原多项式组的解集等于得到的多项式组解集之并
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
# 例: list1={1,2}, list2={a,b}, 则输出 {{1,a},{1,b},{2,a},{2,b}}
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:

# 同上, 不过llist1 可以是 Ltl 输出的多重集合
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:

# 除去 list 中不含未知量的多项式, 除去的因子放在 remfac1 中
Nonumber := proc(list,ord)
	    global remfac1;
            local i,ls;
	    ls := {};
	    for i in list do  
	      if 0 < Class(i,ord) then 
		 ls := ls union {i}
              else remfac1:=remfac1 union {i};
	      fi 
	    od;
	    ls
	    end:

# 列出多项式 pol 的所有不同因子(不考虑重数)
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:

# 输出多项式 p 的初式, 当初式中不含未知量时输出 1
Init := proc(p,ord)
	local pol;
	pol:=Initial(p,ord);
	if Class(pol,ord)=0 then RETURN(1) fi;
	pol
	end:

# 输出多项式组 bset 中的初式的因子(不含重数及不含未知量的因子)
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:

# 输出多项式组 ps 中类等于 1 的多项式个数
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:

# 分解多项式组 rs 中各多项式的因子, 把与 test 中相同的因子加入 remfac 中,
# 输出剩下的因子组(不计因子的重数)
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:

# 把 list 倒过来排
Reverse := proc(list)
	   local i,n,list1;
	   n := nops(list); 
	   list1 := [seq(list[n-i+1],i = 1 .. n)]; 
	   list1
	   end:

# 吴消元法的基本子程序, 假定输入的 ps 不可约的. 输出 {} 表示无解
Cs1 := proc(ps,ord,nzero)
       global AS,iindex;
       local i,j,rm,cset,ps1,bset,rset;
       options remember,system;
       if nums(ps,ord)>1 then RETURN({}) fi;
       rm := {};
       cset := {};
       if {op(ps)} intersect nzero <> {} then RETURN({}) fi;
       bset := basset(ps,ord);
       rset := remseta([op(ps)],bset,ord);
       if rset={1} then Return({}) fi;
       if rset = {} then
	 iindex := iindex+1;
	 AS[iindex] := bset;
	 print(`A New Component`);
	 RETURN({bset})
       fi;
       rset := Nrs(rset,ord,nzero);
       if rset = {} then RETURN({}) fi;
       cset:=map(Cs1,map(`union`,rset,{op(bset)}),ord,nzero);
       map(op,cset);
       end:

Cs2 := proc(ps,ord,nzero)
       local cset,cset1,i,j,rs,rs1;
       options remember,system;
       if nums(ps,ord)>1 then RETURN({}) fi;
       rs1:={};
       cset := Cs1(ps,ord,nzero);
       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) 
       od;
       cset
       end:

Cs := proc(ps,ord,nzero)
      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)});
      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)
      od;
      cset
      end:

Css := proc(ps,ord,nzero)
       global remfac;
       local cset1,cset,i,j,nzero1, ps1,ps2,Is,Ninits;
       options remember,system;
       remfac := {};
       cset1 := Cs(ps,ord,nzero);
       cset := cset1;
       for i in cset1 do
	 Is := Inits(i,ord); 
	 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)
	 od
       od;
       cset
       end:

# 把方程式组 ps 因式分解成分支, 再调用 Css
# remfac1 中含有被除去的常数因子
CSS := proc(ps,ord,nzero)
       global remfac1,iindex;
       local 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;
       remfac1:={};
       iindex:=0;
       cset:={};
       pset:=Nrs(ps,order,nonzero);
       for i to nops(pset) do 
	   cset:=cset union Css(pset[i],order,nonzero,asc) 
       od;
       remfac1:=remfac1 minus nonzero;
       print(`Removed factor(s):`);
       print(remfac1);
       cset;
       end:

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

iindex:=0:
remfac:={}:
remfac1:={};

#---------------------------------------------------
#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:

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:=Ltl_(ls,pss[i])
	    od;
	 fi;
	 ls
	 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]}`,
`       `
):


