wsolve.txt

这是Maple下使用的“吴方法”推理求解多项式方程组的源代码。

end: 

#input  two poly set set1,set2 
#output 
Ltl:=proc(list1,list2) 
local set,i,j; 
     options remember,system; 
   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; 
     options remember,system; 
   set := {}; 
   for i in llist1 do 
       for j in list2 do  set := set union {i union {j}} od 
   od 
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; 
     options remember,system; 
   rss := Produce(rs,ord,test); 
   if rss={} then RETURN({}) fi; 
#    print("rss=");print(rss); 
   rset:=LProduct(rss); 
   rset 
end: 

LProduct:=proc(inlist)    # 输入是集合的集合，输出也是集合的集合 
  local i,j,k,m,n,B,C,D,T; 
     options remember,system;   
   B:=[]; 
   for i from 1 to nops(op(1,inlist)) do 
     B:=[op(B),{op(i,op(1,inlist))}]; 
   end: 
   for i from 2 to nops(inlist) do 
      C:=[];D:=[]; 
      for j from 1 to nops(B) do       
        if ((B[j] intersect op(i,inlist))<>{}) then   
           C:=[op(C),B[j]]; 
         else D:=[op(D),B[j]]; 
        end: 
      end: 
      B:=C; 
      for m from 1 to nops(D) do 
        T:=op(i,inlist); 
        for n from 1 to nops(C) do 
           if (nops(C[n] minus D[m])=1) then 
            T:=T minus (C[n] minus D[m]);   
           end:   
         end: 
       for k from 1 to nops(T) do 
           B:= [op(B),D[m] union {op(k,T)}]; 
        end : 
      end:     
    end:   
    B:={op(B)}; 
 end: 
  
LProduct_wdk:=proc(list1) 
local len,lls,i,j; 
     options remember,system; 
 len:=nops(list1); 
 if len=0 then RETURN(list1) fi; 
 lls:={}; 
 for i in op(1,list1) do 
   lls:=lls union { {i} } 
 od; 
 if len=1 then RETURN (lls) fi; 
 for j from 2 to len do 
   lls:=Lltl(lls, op(j,list1)); 
 od; 
 lls; 
end: 

Lltl:=proc(lls,ls) 
local res,i,j,l1,rm; 
     options remember,system; 
 res:={}; 
 for i in lls do 
   if i intersect ls ={} then 
      for j in ls do 
      res:= res union { i union {j} }; 
      od; 
   else 
      res:=res union {i}; 
   fi; 
 od; 
  
   l1 := nops(res); 
   rm := {}; 
   for i to l1-1 do 
       for j from i+1 to l1 do 
           if op(i,res) minus op(j,res) = {} then 
               rm := rm union {op(j,res)} 
           fi; 
           if op(j,res) minus op(i,res) = {} then 
               rm := rm union {op(i,res)} 
           fi 
       od 
   od; 
   res := res minus rm; 
   res   
end: 

####################################################### 
# Aux functions: 
# 
# 
###################################################### 
Nums:=proc(ps,ord) 
local i,n; 
     options remember,system; 
     n:=0; 
     for i in ps do 
       if Class(i,ord)=1 and POLY_MODE='Algebraic' then n:=n+1 fi ; 
       if Class(i,ord)=1 and (POLY_MODE='Partial_Differential'  or POLY_MODE='Ordinary_Differential') and torder(Leader(i,ord))=0 then n:=n+1 fi; 
     od; 
     n 
end: 

Emptyset:=proc(ds,as,ord) 
local i; 
     options remember,system; 
   for i in ds do if `grobner/normalf`(i,as,ord)=0 then RETURN(1) fi od; 
   0; 
end: 

Non0Constp:=proc(ps,ord) 
local i; 
     options remember,system; 
   for i in ps do if Class(i,ord)=0 and i <> 0 then RETURN(1) fi od; 
   0; 
end:     
    

TestQ:=proc(ps,as,ord) 
global remfac; 
local i; 
     options remember,system; 

   for i in ps do  if `grobner/normalf`(i,as,ord) = 0 then remfac:=remfac union {i};RETURN(1) fi od; 
   0 
end: 

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

JInitials:=proc(bset,ord) 
local pol, product; 
     options remember,system; 
     product:=1; 
     for pol in bset do 
    product:=product*Initial(pol,ord); 
     od; 
     product; 
end: 
        

Inits:=proc(bset,ord) 
local i,list; 
     options remember,system; 
    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: 

#The following will be used in algebraci case ONLY!!! 

Inits1:=proc(bset,ord) 
local i,list; 
     options remember,system; 
    list:={}; 
    for i in bset do if Class(Init(i,ord),ord)>1 then list:=list union Factorlist(Init(i,ord)) fi od; 
###  remove the the factors with class <2 
    for i in list do if Class(i,ord) < 2 then list:=list minus {i} fi od; 
    list 
end: 

####################################################################### 
####################################################################### 
# Compute the Characteristic set with FACTORIZATION      
####################################################################### 
####################################################################### 
# 
# NAME:  Cs_a 
# INPUT: ps:     a polynomial set, Suppose each pol in ps is irreducible over Q. 
#      ord:   indeterminate ordering. if ord:=[z,y,x] means z>y>x; 
#      nzero: a polynomial set. Each pol in nzero is NOT equal to 0 
#                       
#      test:  a polynomial set, which is NOT equal to 0. 
#      T:     a symbol to decide to use which kind of ascending set 
#          T: r_asc, w_asc,g_asc,q_asc=t_asc 
# OUTPUT: a list of ascending set 
##################################################################### 

Cs_a:=proc(ps,ord,nzero,test,T) 
global AS,INDEX,time_bs,time_rs, time_f,__testp__; 
local asc,i,j,rm,cset,test1,ps1,bset,rset, ltime_b,ltime_r,ltime_f; 
options remember,system; 

   if Nums(ps,ord)>1 then RETURN({}) fi; 
   rm := {}; 
   cset := {}; 
  if {op(ps)} intersect nzero <> {} then RETURN({}) fi; 
  if POLY_MODE='Algebraic' then 
       if nops(nzero) <> 0 and  Emptyset(nzero,ps,ord) = 1 then print("An empty set");RETURN({})  fi; 
  
       if __testp__= 1 and nops(test)  <> 0 and TestQ(test ,ps,ord) = 1 then 
           print("One factor of an initial has been found and it will be appended to the original polynomial set "); 
#           print("remfac=",remfac);
           RETURN({}) 
       fi; 
   fi; 
  
#      print("ps=",ps); 
     ltime_b:=time(); 
     bset := Basicset(ps,ord,T); 
#      print("bset=",bset); 
     time_bs:=time_bs+time()-ltime_b; 
     ltime_r:=time(); 
               
     
     rset := Remset([op(ps)],bset,ord,T); 
#       print("rset=",rset); 
     time_rs:=time_rs+time()-ltime_r; 

  
   if rset={1} then RETURN({}) fi; 
# for i in rset do if Class(i,ord)=0 then RETURN({}) fi od; 
   if rset = {} then 
       INDEX:= INDEX+1; 
       AS[INDEX] := bset; 
       print(`A New Component`); 
       RETURN({bset}) 
   fi; 
#    for i in rset do 
#        if (Class(i,ord) = 0) and (i <> 0) then RETURN({}) fi 
#    od; 

############## 
    if POLY_MODE='Algebraic' and __testp__= 1 then
       test1 := test union Inits(bset,ord); 
    else
       test1:=test;
    fi;
##############   
   
   
   ltime_f:=time(); 
   rset := Nrs(rset,ord,nzero union test1); 
   time_f:=time_f+time()-ltime_f; 
   if rset = {} then RETURN({}) fi; 
   cset:=map(Cs_a,map(`union`,rset,{op(bset)}),ord,nzero,test1,T); 
   map(op,cset); 
end: 

####################################################################### 
####################################################################### 
# 
# NAME:  charset_a 
# INPUT: ps:     a polynomial set, Suppose each pol in ps is irreducible over Q. 
#      ord:   indeterminate ordering. if ord:=[z,y,x] means z>y>x; 
#      nzero: a polynomial set. Each pol in nzero is NOT equal to 0 
#                       
#      T:     a symbol to decide to use which kind of ascending set 
#          T: r_asc, w_asc,g_asc,q_asc=t_asc 
# OUTPUT: a list of ascending set 
##################################################################### 

charset_a:=proc(ps,ord,nzero,T) 
global AS,INDEX,time_bs,time_rs, time_f; 
local asc,i,j,rm,cset,test1,ps1,bset,rset, ltime_b,ltime_r,ltime_f; 
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 and Emptyset(nzero,ps,ord) = 1 then print("An empty set");RETURN({})   fi;   

  
#      print("ps=",ps); 
     ltime_b:=time(); 
     bset := Basicset(ps,ord,T); 
#      print("bset=",bset); 
     time_bs:=time_bs+time()-ltime_b; 
     ltime_r:=time(); 
     rset := Remset([op(ps)],bset,ord,T); 
     rset := map(primpart,rset,ord);
#      print("rset=",rset); 
     time_rs:=time_rs+time()-ltime_r; 

  
   if rset={1} or Non0Constp(rset,ord)=1 then RETURN([]) fi; 
# for i in rset do if Class(i,ord)=0 then RETURN({}) fi od; 
   if rset = {} then 
       INDEX:= INDEX+1; 
       AS[INDEX] := bset; 
       print(`A New Component`); 
       RETURN(bset) 
   fi; 
#    for i in rset do 
#        if (Class(i,ord) = 0) and (i <> 0) then RETURN({}) fi 
#    od; 
#    test1 := test union Inits(bset,ord); 
   ltime_f:=time(); 
#    rset := Nrs(rset,ord,nzero union test1); 
   time_f:=time_f+time()-ltime_f; 
#    if rset = {} then RETURN({}) fi; 

   cset:=charset_a({op(bset)} union rset,ord,nzero,T); 
   cset; 
end: 

charset_b:=proc(ps,ord,nzero,T) 
local cset,rs; 
options remember,system; 

   cset := charset_a(ps,ord,nzero,T); 
   rs:=Remset([op(ps)],cset,ord,T);
   rs:=map(primpart,rs,ord);

   if rs={} then RETURN(cset) fi; 
   if rs={1} then RETURN([]) fi; 
   if cset=[] then RETURN([]) fi; 
#   lprint("kkkkkkkkkkkkk");
#   cset:=charset_b({op(ps)} union {op(cset)} union rs,ord,nzero,T); 
   while rs<> {} do
      cset:=charset_a({op(cset)} union rs, ord,nzero,T);
#####we should consider the the following special case which  cset=[]      
      if nops(cset)=0 then RETURN(cset) fi;
      rs:=Remset(ps,cset,ord,T);
   od;
   cset 
end: 

CS_a:=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(Cs_a,pset,order,nonzero,{},asc); 
   [op(map(op,cset))]; 
end: 

####################################################################### 
# 
# NAME:  Cs_b 
# INPUT: ps:     a polynomial set, Suppose each pol in ps is irreducible over Q. 
#      ord:   indeterminate ordering. if ord:=[z,y,x] means z>y>x; 
#      nzero: a polynomial set. Each pol in nzero is NOT equal to 0 
#                       
#      test:  a polynomial set, which is NOT equal to 0. 
#      T:     a symbol to decide to use which kind of ascending set 
#          T: r_asc, w_asc,g_asc,q_asc=t_asc 
# OUTPUT: a list of ascending set 
##################################################################### 

Cs_b:=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 := Cs_a(ps,ord,nzero,test,T); 
   cset1:=cset; 
   for i in cset1 do rs:=Remset([op(ps)],i,ord,'std_asc'); 
   
#    for i in cset1 do rs:=Remset([op(ps)],i,ord,T); 
#                    if rs<>{} then if rs <> {1} then rs1:=rs1 union {rs union {op(i)} } fi; 
if rs<>{} then if rs <> {1} then rs:=Nrs(rs,ord,nzero union test); rs1:=rs1 union map(`union`,rs , {op(i)} ) fi; 
                                   cset:=cset minus {i} fi 
                  od; 
   if rs1={} then RETURN(cset) fi; 
   for j in rs1 do  cset:=cset union Cs_b(ps union j,ord,nzero,test,T) od; 
   cset 
end: 

####################################################################### 
# 
# NAME:  Cs_c 
# INPUT: ps:     a polynomial set, Suppose each pol in ps is irreducible over Q. 
#      ord:   indeterminate ordering. if ord:=[z,y,x] means z>y>x; 
#      nzero: a polynomial set. Each pol in nzero is NOT equal to 0 
#                       
#      test:  a polynomial set, which is NOT equal to 0. 
#      T:     a symbol to decide to use which kind of ascending set 
#          T: r_asc, w_asc,g_asc,q_asc=t_asc 
# OUTPUT: a list of ascending set 
##################################################################### 
Cs_c:=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:=Cs_b({op(ps)},ord,{op(nzero)},{},T); 
    remf:=remfac minus {op(nzero)}; 
#    print(remf);
#     print(remf); 
  for i to nops(remf) do 
#     print(remf[i]); 
#     print(ps); 
       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_c(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:=realfac({op(nzero)},ord) fi; 
   if nargs < 4 then asc:='std_asc' else asc:=T fi; 
   cset:={}; 
   pset:=Nrs(ps,order,nonzero); 
   for i to nops(pset) do cset:=cset union Cs_c(pset[i],order,nonzero,asc) od; 
   [op(cset)] 
end: 

checknum:=0: 
Css:=proc(ps,ord,nzero,T) 
global remfac,checknum; 
local cset1,cset,i,j,nzero1, ps1,ps2,Is,Ninits; 
options remember,system; 
    checknum:=checknum+1; 
   remfac := {}; 
   cset1 := Cs_c(ps,ord,nzero,T); 
   cset := cset1; 
#          nn:=0;
      for i in cset1 do 
#      	  print(nops(cset1));
#      	  nn:=nn+1;