wsolve.txt

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

######################################################### 
#11111111111111111111111111111111111111111111111111111111 
# 
#Part I: The basic functions 
#        Class, Leader, Initial 
######################################################### 

readlib(factors): 
with(Groebner):

      
# the class of poly f wrt variable ordering ord 
Class := proc(dp,ord) 
     local i,idt; 
     options remember,system; 
                     idt:=dIndets(dp); 
                     for i to nops(ord) do 
                         if member(ord[i], idt) then RETURN(nops(ord)-i+1) fi 
                     od; 
                     0 
                 end: 
# the highest order of the indeterminate var in dp 
# INPUT: 
#     dp : a diff pol 
#     var: an indeterminate 
# OUTPUT: the order 
# 
dOrder:=proc(dp,var) 
    local idt,i,dvar; 
    options remember,system;     
    idt:=indets(dp); 
    dvar:=-1;      
    for i in idt do 
      if has(i,var) then 
           if  dvar=-1 or type(dvar,symbol) 
               or (   ( not type(i,symbol) ) and ( op(1,i) > op(1,dvar) ) ) then 
               dvar:=i  ; 
           fi 
     fi; 
    od; 
    if dvar=-1 then -1 
       elif   type(dvar,symbol) then 0 
       else if nops(dvar) =1 then op(1,dvar) 
                else [op(1..nops(dvar),dvar)] 
             fi; 
                 fi; 
    end:     
       
                         
dIndets:=proc(dp) 
    local i,indet_d,idt; 
    options remember,system;      
    idt:={}; 
    indet_d:=indets(dp); 
     
    if POLY_MODE='Algebraic' then RETURN(indet_d) fi; 
     
    for i in indet_d do 
      if type(i,symbol) then idt:={op(idt),i} else idt:={op(idt),op(0,i)} fi; 
    od; 
      idt 
    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,Leader(f,ord)); 
       fi 
   end: 

# the separant of poly f wrt ord 
Separant := 

   proc(f,ord) 
   options remember,system; 
       if Class(f,ord) = 0 then 1 
       else diff(f,Leader(f,ord)); 
       fi 
   end: 

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

# the set of all nonconstant factors of separants of polys in as 
    
Septs:=proc(bset,ord) 
local i,list;
  options remember,system; 
    list:={}; 
    for i in bset do list:=list union Factorlist(Sept(i,ord)) od; 
    for i in list do if Class(i,ord)=0 then list:=list minus {i} fi od; 
    list 
end: 
  
    


Lessp:= 

proc(f,g,ord) 
local cf,cg,cmp,df,dg,leadf,leadg; 
options remember,system; 
   cf := Class(f,ord); 
   cg := Class(g,ord); 
   if type(f,integer) then 
         true 
   elif type(g,integer) then 
        false      
   elif cf = 0 then 
        true     
   elif cg = 0 then 
        false     
   else 
        leadf:=Leader(f,ord); 
        leadg:=Leader(g,ord); 
        
        cmp:=dercmp(leadf,leadg,ord); 
        if cmp = 1 then 
            true 
        elif cmp=-1 then 
            false 
        else         
            df := degree(f,leadf); 
            dg := degree(g,leadg); 
            if df < dg then 
               true 
            elif df = dg then 
               Lessp(Initial(f,ord),Initial(g,ord),ord) 
            else 
               false 
            fi 
        fi 
   fi 
end:     

Least:= 

proc(ps,ord) 
local i,lpol; 
options remember,system; 

  lpol:=op(1,ps); 
  for i from 2 to nops( ps ) do     
    if Lessp(op(i,ps),lpol,ord) then lpol:=op(i,ps) fi; 
  od; 
  lpol 
end: 


#############################################4 Reduced Definitions########################3 

Reducedp:= 
proc(p,q,ord,T) 
local mvar,var,idt; 
options remember,system; 

    mvar:=Leader(q,ord); 
    if type(mvar,symbol) then var:=mvar else var:=op(0,mvar) fi; 
     
     if POLY_MODE='Algebraic' then 
###############Algebraic mode#################################       
       if   nargs <4 or T='std_asc' then 
      if Class(p,ord) > Class(q,ord) and degree(p,mvar) < degree(q,mvar) then 1 else 0 fi; 
    elif T='wu_asc' then 
      if Class(p,ord) > Class(q,ord) and degree(Initial(p,ord),mvar) < degree(q,mvar) then 1 else 0 fi;        
    elif T='gao_asc' then 
    ###########################in the following case , q is a list################# 
      if Class(p,ord) > Class(op(1,q),ord) and Premas(Initial(p,ord),q,ord,'std_asc') <>0 then 1 else 0 fi;      
    elif T='tri_asc' then 
      if Class(p,ord) > Class(q,ord) then 1 else 0 fi;      
    fi;             
###############Algebraic mode#################################          
     elif POLY_MODE='Ordinary_Differential' then 
###############ODE mode#################################           
       if nargs <4 or T='std_asc' then 
        
           if Class(p,ord) > Class(q,ord) and degree(p,mvar) < degree(q,mvar) 
              and dOrder(p,var) <= dOrder(q,var) then 
              1 
           else 
              0 
           fi; 
    elif T='wu_asc' then 
      if Class(p,ord) > Class(q,ord) and degree(Initial(p,ord),mvar) < degree(q,mvar) then 1 else 0 fi;        
    elif T='gao_asc' then 
    ###########################in the following case , q is a list################# 
      if Class(p,ord) > Class(op(1,q),ord) and Premas(Initial(p,ord),q,ord,'std_asc') <>0 then 1 else 0 fi;      
    elif T='tri_asc' then 
      if Class(p,ord) > Class(q,ord) then 1 else 0 fi;                  
            
       fi; 
###############ODE mode#################################     
     elif POLY_MODE='Partial_Differential' then         
###############PDE mode#################################       
       if nargs <4 or T='std_asc' then 
              idt:=indets(p); 
              
           if dercmp(Leader(p,ord), mvar, ord)=-1 and degree(p,mvar) < degree(q,mvar) 
              and   (not member(true, map(Is_proper_derivative,idt,mvar)) )  then 
              1 
           else 
              0 
           fi; 
    elif T='wu_asc' then 
      if Class(p,ord) > Class(q,ord) and degree(Initial(p,ord),mvar) < degree(q,mvar) then 1 else 0 fi;        
    elif T='gao_asc' then 
    ###########################in the following case , q is a list################# 
      if Class(p,ord) > Class(op(1,q),ord) and Premas(Initial(p,ord),q,ord,'std_asc') <>0 then 1 else 0 fi;      
    elif T='tri_asc' then 
      if Class(p,ord) > Class(q,ord) then 1 else 0 fi;                  
            
       fi; 

      
###############PDE mode#################################       
     fi; 
end:      

############################# 
# 
# Name:          Reducedset 
# Input:     ps:     a poly set 
#          p:     a polynomial 
# OUTPUT:     redset:={q | q is in ps, q is rReduced to p } 
############################# 

Reducedset:= 

proc(ps,p,ord,T) 
local i, redset;   
options remember,system; 
    redset:={}; 
    for i in ps do 
      if Reducedp(i,p,ord,T)=1 then redset:={i,op(redset)} fi; 
    od; 
redset 
end:      

############################################################################## 
############################################################################## 
# the basic set of polyset ps 
Basicset:= 

proc(ps,ord,T) 
local qs,b,bs; 
options remember,system; 
         if nops(ps) < 2 then [op(ps)] 
         else 
             qs:=ps;             
             bs:=[]; 
             while (nops(qs) <>0) do 
                b:=Least(qs,ord); 
                bs:=[b,op(bs)]; 
                if T='gao_asc' then 
                  qs :=Reducedset(qs,bs,ord,T);                 
                else 
                  qs :=Reducedset(qs,b,ord,T); 
                fi; 
             od; 
             bs                               
         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,gcd0; 
   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 
#                gcd0:=gcd(l,coeff(r,x,dr),'lu','lv'); 
         gcd0:=gcd(l,coeff(r,x,dr)); 
               lu:=simplify(l/gcd0); 
               lv:=simplify(coeff(r,x,dr)/gcd0); 
         lu:=l; 
         lv:=coeff(r,x,dr); 
               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,T) 
   local r0,r1,asc;
   options remember,system;
   if nargs < 4 then asc:='std_asc' else asc:=T fi;      
     r0:=Premas_a(f,as,ord,asc);
   if POLY_MODE <> 'Partial_Differential' then RETURN(r0) fi;
   
     r1:=Premas_a(r0,as,ord,asc);
     while(r1 <> r0) do
        r0:=r1;
        r1:=Premas_a(r0,as,ord,asc);
     od;
     r1:
end:     
     




Premas_a := 

   proc(f,as,ord,T) 
   local remd,i,degenerate;
   options remember,system;     
       remd := f; 
        
       if nargs <4 then 
        
         for i to nops(as)  do 

            
           remd := NewPrem(remd,as[i],Leader(as[i],ord)); 

            
         od; 
          
       elif T='std_asc' then 
  
         for i to nops(as)  do 
        
           remd := NewPrem(remd,as[i],Leader(as[i],ord)); 
        
         od; 
          
       elif T='tri_asc' then 
        
         for i to nops(as) do 
               if Class(remd,ord) = Class(as[i],ord) then                   
                       remd := NewPrem(remd,as[i],Leader(as[i],ord)); 
               fi 
         od;         
         

###########for wu ascending set#############         
       elif T='wu_asc' then 
        
         for i to nops(as) do 
             remd := WuPrem(remd,as[i],Leader(as[i],ord),ord,'degenerate'); 
             if degenerate=1 then RETURN( Premas(remd,as,ord,T)) fi;
         od;       
##################################################                  
        
       elif T='gao_asc' then 
        
         if Premas( Initial(remd,ord), as,ord,'std_asc' ) =0 then RETURN( Premas(remd,as,ord,'std_asc') ) fi;         
            for i to nops(as) do 
               if Class(remd,ord) = Class(as[i],ord) then 
                       remd := NewPrem(remd,as[i],Leader(as[i],ord));
                       if Premas( Initial(remd,ord), as,ord,'std_asc' ) =0 then RETURN( Premas(remd,as,ord,'std_asc') ) fi;  
 
               fi 
            od;         
          
       fi; 
        
        
       if remd <> 0 then numer(remd/lcoeff(remd)) else 0 fi 
   end: 

    
################################################################################### 
################################################################################### 

# set of non-zero remainders of polys in ps wrt ascending set as 
Remset := 

  proc(ps,as,ord,T) 
  local i,set,pset,r;
  options remember,system; 
  
     set:={}; 
     pset:={op(ps)} minus {op(as)}; 
      
     if POLY_MODE='Partial_Differential' then pset := pset union spolyset(as,ord) fi; 
      
#       print("pset=",pset); 
      
      for i in pset do r:=Premas(i,as,ord,T); 
            if r <>0  and Class(r,ord) =0 
                 then RETURN({1}) 
                 else set:=set union {r} fi 
                  od; 
    set minus {0} 
    end: 
    

################################################################################### 
# Factor the polynomials in Remaider set 
################################################################################### 

Produce:=proc(rs,ord,test) 
global remfac;
options remember,system;  
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: 

Nonumber:=proc(list,ord) 
local i,ls;
  options remember,system;  
   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; 
     options remember,system; 
   list1 := factors(pol)[2]; 
   list2 := {}; 
   for i in list1 do  list2 := list2 union {numer(i[1])} od; 
   list2