wsolve.txt

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

     
    dp:=p; 
    dq:=q; 
     
    for i to nops(ml) do 
      dp:=dDiff(dp,Diff_Var[i],ml[i]-l1[i]); 
      dq:=dDiff(dq,Diff_Var[i],ml[i]-l2[i]); 
    od; 
     
######直觉觉得可以直接用PDPrem(dp,q,Leader(q,ord)) 更好。而这里用的是标准的方法。      
   NPremO(dp,dq, Leader(dq,ord)); 
    
end:     
     
     
maxlist:=proc(l1,l2)      
local i,list;
options remember,system;  
    list:=[]; 
    for i to nops(l1) do 
      list:=[op(list), max(l1[i],l2[i])]; 
    od: 
list 
end:      

#torder   Is_deriv Reducedp Leader Head Leader PDE map Num spoly NewPrem 

# the class of poly f wrt variable ordering ord 
MainVariables := proc(ps,ord) 
    local i,set; 
       options remember,system; 
    set:={}; 
    for i in ps do 
      set:=set union {Leader(i,ord)}; 
    od; 
set 
end:      
            
PRegularp:=proc(p,as,ord) 
    local i,idts,mvs,initp;
    options remember,system;      
    initp:=Initial(p,ord); 
    idts:=indets(initp); 
    if idts={} then RETURN(true) fi; 
    mvs:=MainVariables(as,ord); 
    for i in idts do 
       if member(i,mvs) then RETURN(false) fi; 
    od; 
true 
end:      

ASRegularp:=proc(ps,ord) 
    local i,l,as,p; 
    options remember,system; 
    l:=nops(ps); 
    if l=1 then RETURN(true) fi; 
    as:=[ps[l]]; 
    for i from l-1 to 1 by -1 do 
      p:=ps[i]; 
      if not PRegularp(p,as,ord) then RETURN(false) fi; 
      as:=[p,op(as)]; 
    od;        
true 
end:      

Dec:=proc(ps,ord,nzero,T) 
global remfac,checknum; 
local i,cset1,cset;
options remember,system; 

    cset1:=Cs_c(ps,ord,nzero,T); 
    cset:={}; 
    for i in cset1 do 
      if ASRegularp(i,ord) then 
         cset:=cset union RegDec1(ps,i,ord,nzero,T);      
         else 
                cset:=cset union RegDec2(ps,i,ord,nzero,T); 
         fi;                     
   od; 
cset 
end: 

RegDec1:=proc(ps,as,ord,nzero,T) 
local     cset,j,nzero1, ps1,ps2,Is,Ninits;
options remember,system;

        cset:= {as};      
         if Class(as[nops(as)],ord)=1 then Is:=Inits1(as,ord) else Is := Inits(as,ord) fi;             
         Ninits := Is minus {op(nzero)}; 
     
         for j to nops(Ninits) do 
               ps1 := ({op(ps)} union {op(as)}) union {Ninits[j]}; 
               if j = 1 then nzero1 := {op(nzero)} 
               else nzero1 := nzero1 union {Ninits[j-1]} 
               fi; 
               cset := cset union Dec(ps1,ord,nzero1,T) 
         od; 
cset 
end:           
          
RegDec2:=proc(ps,as,ord,nzero,T) 
local cset,i,j,l,regas,p,initp,lp,pol,r1,r2,Is,Ninits,l2,ps1,nzero1,mvar,vars,mv,d,f,g;
options remember,system; 
    cset:={}; 
     
    l:=nops(as); 
    regas:=[as[l]]; 
     
    for i from 2 to l do 
      p:=as[l-i+1]; 
      mvar:=Leader(p,ord); 
      d:=degree(p,mvar); 
      if PRegularp(p,regas,ord) then 
         regas:= [p,op(regas)] 
      else 
         initp:=Initial(p,ord); 
         vars:=indets(initp); 
         lp:=nops(regas); 
         ###1 开始111 
             for j to lp do 
                 pol:=regas[j]; 
                 mv:=Leader(pol,ord); 
                 if member(mv,vars) then 
#                  注意Resultant函数 

                 r1:=NResultantE(initp,pol, Leader(pol,ord),'f','g'); 
                 r2:=Premas(r1,regas,ord,'std_asc'); 
                 if r2=0 then 
       
                          if Class(as[nops(regas)],ord)=1 then Is:=Inits1(regas,ord) else Is := Inits(regas,ord) fi;             
                          Ninits := Is minus {op(nzero)}; 
               
                       l2:=nops(Ninits); 
                       nzero1:={op(nzero)}; 

                          for j to l2 do 
                              ps1 := ({op(ps)} union {op(as)}) union {Ninits[j]}; 
                              if j <> 1 then 
                               nzero1 := nzero1 union {Ninits[j-1]} 
                              fi; 
                              cset := cset union Dec(ps1,ord,nzero1,T) 
                          od; 
                          print("kkkkkkkkkk"); 
                          if l2 <> 0 then nzero1:=nzero1 union {Ninits[l2]} fi; 
#                   f;g; 
#                   print("cset=",cset);
#                   print("ps=",ps);
#                   print("as=", as);
#                   print("regas=",regas);
                   
                          cset:=cset union Dec(ps union {op(as)} union {op(regas)} union {f},ord,nzero1,T) union Dec(ps union {op(as)} union {op(regas)} union {initp},ord,nzero1,T); 
                   RETURN(cset);            
                 else 
             
#                  p 变形成新的p 
              p:=expand(f*p+g*pol*mvar^d); 
              p:=Premas(p,regas,ord,'std_asc'); 
                 initp:=Initial(p,ord); 
                 vars:=indets(initp); 
              fi; 
              fi; 
             od;   
         ###1 结束111   
         regas:=[p,op(regas)]   
      fi;        
      od; 
    cset:=RegDec1(ps,regas,ord,nzero,T); 
cset      
end:        


#如果首项系数变成0，而余式不是0则置degenerate为1否则为0。
WuPrem:=

   proc(uu,vv,x,ord,degenerate)
   local r,init,lead,lmono,red,init_r,s,t;
   options remember,system;    
   degenerate:=0;
   if uu=0 then RETURN(0) fi;
      r:=uu;  
   if Class(uu,ord) = Class(vv,ord) then 
      r:= NPrem(uu,vv,x);
   elif Class(uu,ord) > Class(vv,ord) then
      init:=Initial(uu,ord); 
      if  degree(init,x) >0 and Reducedp(init, vv,ord,'std_asc') = 0  then 
          lead:=Leader(uu,ord);
#	  collect(uu,lead);
          lmono:=lead^degree(uu,lead);
          red:=expand(uu-init*lmono);
          init_r:=NPremE(init,vv,x,'s','t');    
          r:=expand(init_r*lmono+s*red); 
          if init_r=0 and r <>0 then degenerate:=1 fi;
#          collect(r,lead);
       fi;
    fi;
    r;
end:    

NPremE := 

   proc(uu,vv,x,s1,t1) 
   local r,v,dr,dv,l,t,lu,lv,s0,t0; 
   options remember,system; 

        s0:=1; 
        t0:=0;           
           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'); 
         s0:=lu*s0; 
         t0:=lu*t0+lv*x^(dr-dv); 
               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; 
        s1:=expand(s0); 
        t1:=expand(t0);           
           r 
   end: 

#This function is for ascending set normalization. 
#in fact, it is not  a resultant for 2 pols with indeterminate x. 
#for two polys uu,vv with indeterminate x 
#it has the property m*uu+n*vv=r with degree(r,x)=0; 
NResultantE := 

   proc(uu,vv,x,m,n) 
   local r,s,t,r0,s0,t0,r1,s1,t1,tmpr,tmps,tmpt; 
   options remember,system; 
        r0:=uu; 
        s0:=1; 
        t0:=0; 
        r1:=vv; 
        s1:=0; 
        t1:=1;           
        if degree(r0,x)=0 then m:=1;n:=0;RETURN(r0) fi; 
           if degree(r1,x)=0 then m:=0;n:=1;RETURN(r1) fi; 

           while degree(r1,x) >= 1 do 
             r:=NPremE(r0,r1,x,'s','t');        
          
            tmpr:=r1; 
         tmps:=s1; 
         tmpt:=t1; 
         r1:=r; 
         s1:=s*s0-t*s1; 
         t1:=s*t0-t*t1; 
         r0:=tmpr; 
         s0:=tmps; 
         t0:=tmpt; 
           od; 
        m:=s1; 
        n:=t1; 
           r 
   end: 

#printf("    Pls type 'with(linalg)' firstly before calling function 'Resultant_E' "); 

##The following is to compute the resultant of 2 pols via computing the determinant directly.
# input: M is a list of lists , whereby to denote a matrix; 
#        i,j are two integer numbers; 
# output:sM is a new list of lists, obtained by removing 
#        the ith row and jth column of L, 
SubM:=proc(M,i,j) 
    local row,tmp,sM; 
    if i<1 or i>nops(M) 
         then error "Input row index %1 broke bounds",i;fi; 
    if j<1 or j>nops(M[1]) 
         then error "Input column index %1 broke bounds",j;fi; 
    sM:=[]; 
    tmp:=subsop(i=NULL,M); 
    for row in tmp 
         do row:=subsop(j=NULL,row); 
            sM:=[op(sM),row]; 
         od; 
    sM; 
end: 
     



# with(linalg) firstly 
# input: A,B are two polynomials, v is a variable; 
# output:[R,T,U]; R,T,U are three polynomials, where R is the resultant 
#        of A,B w.r.t v, so that A*T+B*U=R 
#with(linalg); 
NResultantE_Matrix:=proc(A,B,v,TA,UB) 
    local m,n,d1,d2,cA,cB,row,i,j,syl,msyl,sM,sign,R,T,U,result; 
#    lprint("nargs=",nargs); 

    if (nargs <> 3) and (nargs <> 5) then error "Number of parameters should be 3 or 5" fi; 
     
    m:=degree(A,v); 
    n:=degree(B,v); 
    if m=0 then error "InPut Polynomial %1 Has No Variable %2",A,v;fi; 
    if n=0 then error "InPut Polynomial %1 Has No Variable %2",B,v;fi; 
  
#  to get coefficients of A and B 
    cA:=[]; 
    for d1 from 0 to m 
         do cA:=[coeff(A,v,d1),op(cA)]; 
         od; 
    cB:=[]; 
    for d2 from 0 to n 
         do cB:=[coeff(B,v,d2),op(cB)]; 
         od; 

#  initialize sylvester matrix 
     
    syl:=[]; 
    for i from 1 to n         
         do row:=cA; 
               for j from 1 to i-1 
              do row:=[0,op(row)];od; 
            for j from 1 to n-i 
                do row:=[op(row),0];od; 
               syl:=[op(syl),row]; 
            od; 
             
    for i from 1 to m 
         do row:=cB; 
               for j from 1 to i-1 
              do row:=[0,op(row)];od; 
            for j from 1 to m-i 
                do row:=[op(row),0];od; 
               syl:=[op(syl),row]; 
            od; 
     
    msyl:=linalg[matrix](m+n,m+n,syl); 
     
#   to get resultant of A,B w.r.t v 
    R:=linalg[det](msyl); 
    R:=expand(R); 
    if nargs=3 then RETURN(R) fi; 
     
#   to get T 
    T:=0; 
    for d1 from 1 to n 
         do 
         sign:=(-1)^(d1+m+n mod 2);           
         sM:=SubM(syl,d1,m+n); 
         T:=T+(v^(n-d1))*expand(sign*linalg[det](linalg[matrix](m+n-1,m+n-1,sM))); 
         T:=expand(T); 
         od; 
          
#  to get U 
        U:=0; 
    for d2 from 1 to m 
         do 
         sign:=(-1)^(d2+m mod 2);      
         sM:=SubM(syl,d2+n,m+n); 
         U:=U+(v^(m-d2))*expand(sign*linalg[det](linalg[matrix](m+n-1,m+n-1,sM))); 
         U:=expand(U); 
         od; 
        if nargs=5 then TA:=T; UB:=U fi; 
    R; 
end: 


############## computing partioned-parametric Grobner basis###################
interface(warnlevel=0);
with(Groebner): 


########p在约束Zero(a11/a12)下一定不为０,ord是参变元########## 
NotZeroUnderCons:=proc(a11,a12,p,ord) 
option remember; 
    RETURN(IsEmpty({op(a11),p},a12,ord,ord)); 
end: 

########p在约束Zero(a11/a12)下一定为０,ord是参变元########## 
IsZeroUnderCons:=proc(a11,a12,p,ord) 
option remember; 
    RETURN(IsEmpty(a11,{op(a12),p},ord,ord)); 
end: 

########约束多项式conspolyset形式:[[{a11}，{a12}]，[{a21}，{a22}]]##
UnAmbiguous:=proc(conspolyset,var_para::list,term_order) 
local p,vp,result,a11,a12,a21,a22,newa22,lc,vars,paras,reduct; 

vars:=var_para[1]; 
paras:=var_para[2]; 

a11:=conspolyset[1][1]; 
a12:=conspolyset[1][2]; 
a21:=conspolyset[2][1]; 
a22:=conspolyset[2][2]; 
if a22={} then RETURN ({conspolyset}); fi; 

p:=a22[1]; 
newa22:=a22 minus {p};

p:=numer(normalf(p,a21,term_order(op(vars))));

if nops(a11) <>0 then p:=numer(normalf(p,a11,term_order(op(paras)))); fi;

#########如果p是0多项式############ 
if p=0 then  result:=UnAmbiguous([[a11,a12],[a21,newa22]],var_para,term_order);  RETURN(result); fi; 

##########如果p是非0常数##############   
if (p <> 0) and ( type(p,'constant') ) then RETURN ( {[conspolyset[1],[{1},{} ]  ]}); fi; 
  
lc:=leadcoeff(p,term_order(op(vars))); 
reduct:=expand(p-lc*leadterm(p,term_order(op(vars)))); 
vp:=indets(p); 
    
#######如果p含有变元############   
 if (vp intersect {op(vars)}) <> {} then 
     if(type(lc,'constant')=true) then 
        result:=UnAmbiguous([[a11,a12],[a21 union {p}, newa22]],var_para,term_order); 

     elif NotZeroUnderCons(a11,a12,lc,paras) 
         then result:=UnAmbiguous([[a11,a12],[a21 union {p}, newa22]],var_para,term_order); 
             
     elif IsZeroUnderCons(a11,a12,lc,paras) 
         then result:=UnAmbiguous([[a11,a12],[a21,newa22 union {reduct}]],var_para,term_order); 
           
     else         
         lc:=CutNozeroFactors(lc,a12);            
         result:=UnAmbiguous([[{op(a11),lc},a12],[a21,newa22 union {reduct}]],var_para,term_order) 
               union UnAmbiguous([[a11,a12 union FactorList(lc)],[a21 union {p}, newa22]],var_para,term_order); 
                 
     fi; 
 else 

#######如果p不含变元，只含参数############   
##modify 
     if NotZeroUnderCons(a11,a12,p,paras) then result:={[conspolyset[1],[{1},{} ]]}; 
      
     elif IsZeroUnderCons(a11,a12,p,paras)  then result:=UnAmbiguous([[a11,a12],[a21,newa22]],var_para,term_order) ;           
        
     else 
          lc:=CutNozeroFactors(lc,a12);
          result:=UnAmbiguous([[{op(a11),lc},a12],[a21,newa22]],var_para,term_order) 
                  union UnAmbiguous([[a11,a12 union FactorList(lc)],[{1}, {}]],var_para,term_order); 
     fi; 
 fi; 

result: 
end: 

######pol对list1中的多项式逐个做除法，返回商#######
DivideList:=proc(pol,list1)
local i,q,rt;
   rt:=pol;
   for i in list1