点击打开链接 王定康研究员的中科院 数学与系统科学研究院 主页;

通常  x:\program files\maple xx\ 的路径是

currentdir() 命令返回的路径, 下载 点击打开链接 文本格式代码文件到"当前路径"即可

# #1 # #Part I: The basic functions # Class, Leader, Initial # readlib(factors): with(Groebner): _field_characteristic_:=0; # 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 Definitions3 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: # 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 NPremfinite := proc(uu,vv,x,ord) local r,v,dr,dv,l,t,lu,lv,gcd0; options remember,system; global _field_characteristic_; 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) mod _field_characteristic_; if dr = 0 then r := 0 else r := subs(x^dr = 0,r) mod _field_characteristic_ fi; r := expand(lu*r)-t mod _field_characteristic_; # r :=finite_simplify(r,ord,_field_characteristic_); dr := degree(r,x) od; r fi end: finite_simplify:=proc(r,ord,p) local i,res; res:=r; for i to nops(ord) do res:=subs(ord[i]^p=ord[i],res) mod p; od; res: 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; global _field_characteristic_ ; remd := f; 下面专门处理有限域的情形 if (_field_characteristic_ <> 0) then for i to nops(as) do remd := NPremfinite(remd,as[i],Leader(as[i],ord),ord); od; fi; 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: misc procedures Reverse:=proc(list) local i,n,list1; n := nops(list); list1 := [seq(list[n-i+1],i = 1 .. n)]; list1 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; if (_field_characteristic_ <>0) then list1:=(Factors(pol) mod _field_characteristic_)[2]; else list1 := factors(pol)[2]; fi; list2 := {}; for i in list1 do list2 := list2 union {numer(i[1])} od; list2 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; if {op(ds)} intersect {op(as)} <> {} then RETURN(1) fi; # for i in ds do if grobner[normalf](i,as,tdeg(op(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 asc_set,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; asc_set[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 asc_set,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; asc_set[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; # print(nn); if Class(i[nops(i)],ord)=1 and POLY_MODE='Algebraic' then Is:=Inits1(i,ord) else Is := Inits(i,ord) fi; if POLY_MODE='Partial_Differential' or POLY_MODE='Ordinary_Differential' then Is := Is union Septs(i,ord) fi; Ninits := Is minus {op(nzero)}; # print("checknum=",checknum, "Ninits=",Ninits); 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: realfac:=proc(ps,ord) local s,res,i; options remember,system; s:=Produce(ps,ord,{}); res:={}; for i in s do res:=res union i; od; res; end: Degenerate:=proc(ds,as,ord) local i,r; options remember,system; for i in ds do r:=Premas(i,as,ord,'std_asc'); if r =0 then return 1 fi; od; 0; end: POLY_MODE:= one of ['Algebraic','Ordinary_Differential','Partial_Differential'] depend_relation is like : [[x,y],[t]]; T:= one of ['std_asc','norm_asc','wu_asc','gao_asc','tri_asc'] CSS:=proc(ps,ord,nzero,T) global factor_time,prem_time, t_time,time_bs,time_rs,time_f; local asc,i,pset,cset,nonzero,order,result,wm; options remember,system; 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; if asc='norm_asc' then return Dec({op(ps)},order,nonzero,'std_asc');fi; cset:={}; factor_time:=0; prem_time:=0; time_bs:=0; time_rs:=0; time_f:=0; t_time:=time(); pset:=Nrs(ps,order,nonzero); if asc='norm_asc' then for i to nops(pset) do cset:=cset union Dec(pset[i],order,nonzero,'std_asc') od; else for i to nops(pset) do cset:=cset union Css(pset[i],order,nonzero,asc) od; fi; result:=[]; for i in cset do if Degenerate(nonzero,i,order)<>1 then result:=[op(result),i]; fi od; #lprint("Timing","Total", time()-t_time, "Factoring", factor_time,"Prem",prem_time); #lprint("BasicSet",time_bs,"RS",time_rs,"factor",time_f); result; end: charset:=proc(ps,ord,nzero,T) global factor_time,prem_time, t_time,time_bs,time_rs,time_f; local asc,nonzero,order,result,cm; options remember,system; 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; prem_time:=0; time_bs:=0; t_time:=time(); result:=charset_b({op(ps)}, order,nonzero,asc); # result:=charset_b(map(primpart,{op(ps)},order),order,nonzero,asc); #lprint("Timing","Total", time()-t_time, "BasicSet",time_bs,"Prem",prem_time); result; end: Css1:=proc(ps,ord,test,type) local cset,i,j,septs,cset1,nonzero,vv; options remember,system; cset:=CSS(ps,ord,test,type); # septs:={}; cset1:=cset; nonzero:={op(test)}; for i in cset1 do septs:=Septs(i,ord) minus Inits(i,ord); if septs<>{} then for j in septs do vv:= {op(Css1(ps union {op(i)} union {j},ord,nonzero,type))}; nonzero:={op(nonzero),j}; cset:={op(cset)} union vv ; od fi od; [op(cset)] end: Gsolve:=proc(ps,ord,test) gsolve(ps,ord,test) end: e_val:=proc(ps,ord,test,T) global factor_time,prem_time, t_time,time_bs,time_rs,time_f; local asc,nonzero,order,result,cm; options remember,system; if nargs < 2 then error("The input should have at least 2 parameters") fi; if nargs < 3 then nonzero:={} else nonzero:=realfac({op(nzero)},ord) fi; if nargs < 4 then asc:='std_asc' else asc:=T fi; prem_time:=0; time_bs:=0; t_time:=time(); nonzero:=test; result:=Css1({op(ps)}, ord,nonzero,asc); #lprint("Timing","Total", time()-t_time, "BasicSet",time_bs,"Prem",prem_time); result; end: wsolve:=proc() local rt; global POLY_MODE; POLY_MODE:='Algebraic'; rt:=CSS(args); rt; end: od_wsolve:=proc() local rt; global POLY_MODE,depend_relation; POLY_MODE:='Ordinary_Differential'; if type(depend_relation,'list') and nops(depend_relation) >1 and type(depend_relation[2],'list') and nops(depend_relation[2])=1 then rt:=CSS(args); else #lprint("Please set depend_relation:=[[x,y,...],[t]] first, if x,y,... depend on t"); rt:=0; fi; rt; end: pd_wsolve:=proc() local rt; global POLY_MODE,depend_relation; POLY_MODE:='Partial_Differential'; if type(depend_relation,'list') and nops(depend_relation) >1 and type(depend_relation[2],'list') and nops(depend_relation[2])>1 then rt:=CSS(args); else #lprint("Please set depend_relation:=[[x1,x2,...],[t1,t2,...]] first, if x1,x2,... depend on t1,t2,..."); rt:=0; fi; rt; end: INDEX:=0: __testp__:=0: remfac:={}: prem_time:=0: factor_time:=0: time_bs:=0: time_rs:=0: time_f:=0: POLY_MODE:=Algebraic: RANKING_MODE:='cls_ord_var': # Test Program for Differentiation #---------------------------------------- # Global dependence #---------------------------------------- # #1 # #Part IV: The basic functions for differential poly # # #Diff_Var:=[u,v,w]; # Declare variables; #Diff_Indets:=[X1,X2,X3]; # Declare Diff ndeterminates; #Lvar:=nops(Diff_Var); #---------------------------------------------------------- # main function #---------------------------------------------------------- dDiff:=proc(dp,var,n) local df,i; options remember,system; df:=dp; if nargs=2 then RETURN (DF(dp,var)) fi; for i to n do df:=DF(df,var); od; df; end: #========================================================== # # dq <- DF(dp, var) # # (differetiating a diff polynomial) # # # Input: dp, a diff polynomial; # var, the variable w.r.t which to be differentiate # # Output: dq, the result after differentiating dp w.r.t var; # #=========================================================== DF:=proc(dp, var) local dq, dp1, dp2; options remember,system; #Step1 [Recursive base] if op(0,dp) <> `+` then dq:=DF_prod(dp,var); RETURN(dq); fi; #Step2 [Recursion] dp1:=op(1,dp); dp2:=normal(dp-dp1); dq:=normal(DF_prod(dp1,var)+DF(dp2,var)); #Step3 [Prepare for return] RETURN(dq); end: #------------------------------------------ # subroutines #------------------------------------------ #========================================================= # # # der <- DF_indet(indet, var) # # (differentiate a differential indeterminante) # # Input: indet, a differential indeterminate which is a # derivative of symobls in Diff_Indets w.r.t # some variables in Diff_Var; # # var, the variable w.r.t which to be differeniated # # Output: der, the derivative of indet w.r.t var # #=========================================================== DF_indet:=proc(indet, var) local der, p, i, index, j,Diff_Var,Diff_Indets,Lvar; global depend_relation; options remember,system; #STEP0 [Initiation diffs]; Diff_Var:=depend_relation[2]; Diff_Indets:=depend_relation[1]; Lvar:=nops(Diff_Var); #STEP1 [Special cases] if not member(var, Diff_Var, 'p') then der := diff(indet, var); RETURN(der); fi; if not member(indet, Diff_Indets) and not member(op(0,indet), Diff_Indets) then der := diff(indet, var); RETURN(der); fi; #STEP2 [Compute] index:=NULL; #Initialization if member(indet, Diff_Indets) then # build a derivative from a diff indet for i from 1 to Lvar do if i = p then index:=index,1; else index:=index,0; fi; od; der:=indet[index]; else # modify a derivative from a given one for i from 1 to Lvar do j:=op(i,indet); if i = p then index:=index,(j+1); else index:=index,j; fi; od; der:=op(0,indet)[index] fi; #STEP3 [Prepare for return] RETURN(der); end: #======================================================== # # dq <- DF_power(dp, var) # # (differentiating a power of a drivative of a diff indet) # # input: dp, a power of a diff indet # var, the variable w.r.t which to be differentiated # # output: dq, the result after differentiating dp w.r.t var # #========================================================= DF_power:=proc(dp, var) local dq, indet, expon; options remember,system; #Step1 [Special cases] # constant case if dp = 1 then dq := 0; RETURN(dq); fi; # indeterminate case if op(0,dp) <> `^` then dq := DF_indet(dp, var); RETURN(dq); fi; #Step2 [Compute] indet:=op(1,dp); expon:=op(2,dp); dq:=expon*indet^(expon-1)*DF_indet(indet,var); #Step3 [Prepare for return] RETURN(dq); end: #========================================================= # # dq <- DF_prod(dp, var) # # (differentiating a product of derivatives of some diff indets) # # input: dp, a product of derivatives of some diff indets # var, the variable w.r.t which to be differentiated # # output: dq, the result after differentiating dp w.r.t var; # #========================================================== DF_prod:=proc(dp, var) local dq, dp1, dp2; options remember,system; #Step1 [Recursive base] if op(0,dp) <> `*` then dq := DF_power(dp, var); RETURN(dq); fi; #Step2 [Recursion] dp1:=op(1,dp); dp2:=normal(dp/dp1); dq:=normal(DF_power(dp1,var)*dp2+dp1*DF_prod(dp2,var)); #Step3 [Prepare for return] RETURN(dq); end: #=============================================================== # # ord <- torder(der) # # computing the order of a derivative of a diff indet # # input: der, a derivative of a diff indet # # output: ord, ord is the order of der. # #================================================================= torder:=proc(der) local ord,i,Diff_Var,Lvar; global depend_relation; options remember,system; #STEP1 [Initialize] if type(der,symbol) then RETURN(0) fi; ord := 0; #STEP2 [Compute] for i from 1 to nops(der) do ord := ord + op(i, der); od; ord; end: #================================================================== # # lead <- Leader(dp, ranking,ord, _cls, _ord) # # Input: dp, a nonzero diff poly; # ranking, a specified ranking; # # Output: the Leader w.r.t. ranking; # # Option: _cls, the class of the lead; # _ord,, the order of the lead; # #=================================================================== Mainvar:=proc(p,ord) options remember,system; Leader(p,ord); end: Leader := proc(dp, ord, _cls, _ord) local lead, cls, order, L, l, der, clsord, i, j, cls1, ord1,Diff_Var,Lvar; global depend_relation; options remember,system; #Step1 [Initialize] Diff_Var:=depend_relation[2]; Lvar:=nops(Diff_Var); L:=indets(dp); l:=nops(L); lead := 1; cls := 0; order := 0; #Step 2 [Algebraic mode] if POLY_MODE='Algebraic' then for i to nops(ord) do if member(ord[i],L) then RETURN(ord[i]) fi od; fi; #Step2 [cls_ord_var case] if RANKING_MODE = cls_ord_var then for i from 1 to l do der := L[i]; cls1:=Class(der,ord); ord1:=torder(der); if cls1 > cls or (cls1 = cls and ord1 > order) then lead := der; cls := cls1; order := ord1 else if cls1 = cls and ord1 = order and ord1 > 0 then for j from 1 to Lvar do if op(j, lead) < op(j, der) then lead := der; cls := cls1; order := ord1; break; else if op(j, lead) > op(j, der) then; break; fi; fi; od; fi; fi; od; fi; #Step3 [ord_cls_var case] if RANKING_MODE = ord_cls_var then for i from 1 to l do der := L[i]; cls1:=Class(der,ord); ord1:=torder(der); if ord1 > order or (ord1 = order and cls1 > cls) then lead := der; cls := cls1; order := ord1 elif ord1 = order and cls1 = cls and ord1 > 0 then for j from 1 to Lvar do if op(j, lead) < op(j, der) then lead := der; cls := cls1; order := ord1; break; elif op(j, lead) > op(j, der) then; break; fi; od; fi; od; fi; #STEP3 [Prepare for return] if nargs > 2 then _cls := cls; fi; if nargs > 3 then _ord := order; fi; RETURN(lead); end: depend:=proc(l1,l2) global depend_relation; depend_relation:=[l1,l2]; 1: end: #============================================================ # # {1,-1,0} <- dercmp(der1, der2) # # input: der1, der2, two derivatives; # # # output: 0 if der1 = der2 # 1 if der1 < der2 # -1 if der1 > der2 # #============================================================ dercmp:=proc(der1, der2,ord) local dp, lead; options remember,system; if der1=der2 then RETURN(0) fi; dp := der1 + der2; lead:=Leader(dp,ord); if lead = der1 then RETURN(-1); else RETURN(1); fi; end: #if der1 is a derivative of der2 then return true else false Is_derivative:=proc(der1, der2) local dt1,dt2,i, l1,l2; options remember,system; if type(der1,symbol) then dt1:=der1 else dt1:=op(0,der1) fi; if type(der2,symbol) then dt2:=der2 else dt2:=op(0,der2) fi; if dt1<>dt2 then RETURN(false) fi; if der1=der2 then RETURN(true) fi; if torder(der2)=0 then RETURN(true) elif torder(der1)=0 then RETURN(false) ; else l1:=dOrder(der1,dt1); l2:=dOrder(der2,dt1); for i to nops(l1) do if op(i,l1) < op(i,l2) then RETURN(false) fi; od; fi; true; end: #if der1 is a proper derivative of der2 then return true else false Is_proper_derivative:=proc(der1, der2) local dt1,dt2,i, l1,l2; options remember,system; if type(der1,symbol) then dt1:=der1 else dt1:=op(0,der1) fi; if type(der2,symbol) then dt2:=der2 else dt2:=op(0,der2) fi; if dt1<>dt2 then RETURN(false) fi; if der1=der2 then RETURN(false) fi; if torder(der2)=0 then RETURN(true) elif torder(der1)=0 then RETURN(false) ; else l1:=dOrder(der1,dt1); l2:=dOrder(der2,dt1); for i to nops(l1) do if op(i,l1) < op(i,l2) then RETURN(false) fi; od; fi; true; end: #return all the derivatives in dp of der proper_derivatives:=proc(dp,der) local i,idt,dets; options remember,system; idt:=indets(dp); dets:={}; for i in idt do if Is_proper_derivative(i,der) then dets:={i,op(dets)}; fi; od; dets; end: #sometimes, Maple's GCD is NOT valid for parameters which has indexes such a[1],a[2]; NPremO := proc(uu,vv,xx) local r,v,dr,dv,lcv,rtv,g,lu,lv; options remember,system; if degree(vv,xx)=0 then RETURN(0) fi; if type(vv/xx,integer) then coeff(uu,xx,0) else r := expand(uu); dr := degree(r,xx); v := expand(vv); dv := degree(v,xx); lcv:=lcoeff(v,xx); # rtv:=sterm(v,xx); rtv:=v-expand(lcoeff(v,xx)*xx^degree(v,xx)); while dv <= dr and r <> 0 do # g:=gcd(lcv,coeff(r,xx,dr)); lu:=lcv; lv:=coeff(r,xx,dr); r:=expand(lu*(r-expand(lcoeff(r,xx)*xx^degree(r,xx)))-lv*rtv*xx^(dr-dv)); dr := degree(r,xx) od; r fi end: 变量y有两种情形.y and y[1]所以要小心. ODPrem:= proc(uu,vv,y) local u,x,dv,du,d,t,var; global depend_relation; options remember,system; t:=depend_relation[2][1]; var:=depend_relation[1]; # if dClass(uu,var)=0 then RETURN(uu) fi; u:=uu; if type(y,symbol) then x:=y else x:=op(0,y) fi; dv:=dOrder(vv,x); du:=dOrder(u,x); d:=du-dv; if d<0 then RETURN(uu) fi; #the following program is for ordinary differential case while d>0 do u:=NPremO(u,dDiff(vv,t,d), x[du]); du:=dOrder(u,x); d:=du-dv; od; # if dv =0 then NPremO(u,vv,x) else NPremO(u,vv,x[dv]) fi; end: NewPrem:= proc(uu,vv,y) local u,x,dv,du,d,t,var; global depend_relation; options remember,system; if POLY_MODE='Algebraic' then NPrem(uu,vv,y); elif POLY_MODE='Ordinary_Differential' then ODPrem(uu,vv,y); elif POLY_MODE='Partial_Differential' then PDPrem(uu,vv,y); fi; end: Headder:=proc(idts) local hder,i,j; options remember,system; hder:=op(1,idts); # if type(hder,symbol) then var:=hder else var:=op(0,hder) fi; for i in idts do if torder(i) > torder(hder) then hder:=i ; elif torder(i)=torder(her) then for j to nops(hder) do if op(j,i) > op(j,hder) then hder:=i;break; elif op(j,i) < op(j,hder) then break; fi; od; fi; od; hder: end: PDPrem:=proc(dp, dq, y) local Diff_Indets,Diff_Var, var,dr, idt,l2,lder,da,l1,i,s; options remember,system; #Step1 [Special case] Diff_Indets:=depend_relation[1]; Diff_Var:=depend_relation[2]; if y=1 then RETURN(0) fi; if type(y,symbol) then var:=y else var:=op(0,y) fi; dr:=dp; idt:=proper_derivatives(dr,y); l2:=dOrder(y,var); #Step2 [Recursive base] while nops(idt)<>0 do lder:=Headder(idt); da := dq; l1:=dOrder(lder,var); for i to nops(l1) do if l2=0 then s:=op(i,l1) else s:=op(i,l1) -op(i,l2) fi; da := dDiff(da, Diff_Var[i],s); od; dr := NPremO(dr, da, lder); idt:=proper_derivatives(dr,y); od; dr:=NPremO(dr,dq,y) ; #Step4 [Prepare for return] dr; end: # # set: A set; cmb:=proc(set0) local p,q, ls,set1; options remember,system; ls:={}; set1:=set0; for p in set0 do set1:=set1 minus {p}; for q in set1 do ls:={op(ls),[p,q]} od; od; ls; end: # modified Nov.6 # spolyset:=proc(as,ord) local res,bs,l,p,cls,i,poly; options remember,system; res:={}; bs:={op(as)}; l:=nops(bs); if l <= 1 then RETURN(res) fi; p:=op(1,bs); cls:=Class(Leader(p,ord),ord); bs:=bs minus {p}; while nops(bs) <>0 do for i in bs do if Class(Leader(i,ord),ord) = cls then poly:=PD_spoly(i,p,ord); res:=res union {poly} fi; od; p:=op(1,bs); cls:=Class(Leader(p,ord),ord); bs:=bs minus {p}; od; res minus {0}; end: PD_spoly:=proc(p,q,ord) local leadp,leadq,var,l1,l2,ml,dp,dq,i,Diff_Var; global depend_relation; options remember,system; leadp:=Leader(p,ord); leadq:=Leader(q,ord); var:=op(0,leadp); Diff_Var:=depend_relation[2]; l1:=dOrder(leadp,var); l2:=dOrder(leadq,var); ml:=maxlist(l1,l2); 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: PNormalp:=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: ASNormalp:=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 PNormalp(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 ASNormalp(i,ord) then cset:=cset union NormDec1(ps,i,ord,nzero,T); else cset:=cset union NormDec2(ps,i,ord,nzero,T); fi; od; cset end: NormDec1:=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: NormDec2:=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 PNormalp(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:=NormDec1(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: # # Name: SylvesterResultant # Calling sequence: SylvesterResultant(f, g, x, u, v) # Input: f, g two polynomials in x with positive degrees # x, a name # Output: the Sylvester resultant of f and g with respect to x # Optional: u, v, such that # u*f+v*g = resultant #SylvesterResultant := proc(f, g, x, u, v) NResultantE := proc(f, g, x, u, v) local df, dg, deg, A, a, b, U, V, i, m, Q, l, e, s, us, vs; #------------------------------------------- # check input #------------------------------------------- (df, dg) := degree(f, x), degree(g, x); if df = 0 or dg = 0 then error "input polynomials should be of positive degrees"; end if; #--------------------------------------------- # initialize data #--------------------------------------------- if df >= dg then (A[1], A[2]) := f, g; (deg[1], deg[2]) := df, dg; else (A[1], A[2]) := g, f; (deg[1], deg[2]) := dg, df; end if; if nargs > 3 then (U[1], U[2]) := 1, 0; if nargs > 4 then (V[1], V[2]) := 0, 1; end if; end if; (a[1], b[1], a[2], b[2]) := (1, 1, lcoeff(A[2],x), lcoeff(A[2], x)^(deg[1]-deg[2])); i := 2; l[i] := deg[i-1]-deg[i]+1; #---------------------------------------------- # compute #---------------------------------------------- while true do i := i + 1; #------------------------------- # form extraneous factor #------------------------------- e[i] := (-1)^l[i-1]*a[i-2]*b[i-2]; #-------------------------------- # compute remainder #-------------------------------- if nargs = 3 then A[i] := evala(Simplify(prem(A[i-2], A[i-1], x)/e[i])); else A[i] := evala(Simplify(prem(A[i-2], A[i-1], x, 'm', 'Q')/e[i])); U[i] := evala(Simplify((m*U[i-2] - Q*U[i-1])/e[i])); if nargs = 5 then V[i] := evala(Simplify((m*V[i-2] - Q*V[i-1])/e[i])); end if; end if; #----------------------------------- # Resultant is zero #----------------------------------- if A[i] = 0 then if nargs > 3 then u := U[i]; if nargs > 4 then v := V[i]; end if; end if; return A[i]; end if; #------------------------------------ # update auxiliary sequences #------------------------------------ deg[i] := degree(A[i], x); l[i] := deg[i-1]-deg[i]+1; a[i] := lcoeff(A[i], x); b[i] := evala(Simplify(a[i]^(l[i]-1)/b[i-1]^(l[i]-2))); #------------------------------------------- # Resultant is nonzero #------------------------------------------- if deg[i] = 0 then if nargs > 3 then s := evala(Simplify(b[i]/a[i])); us := evala(Simplify(U[i]*s)); if nargs > 4 then vs := evala(Simplify(V[i]*s)); end if; if df >= dg then u := us; if nargs > 4 then v := vs; end if; else u := (-1)^(df*dg)*vs; if nargs > 4 then v := (-1)^(df*dg)*us; end if; end if; end if; if df >= dg then return b[i]; else return (-1)^(df*dg)*b[i]; end if; end if; end do; end proc: #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 2005.11.4 需要　“wsolve” 和　“project” 文后付有主要函数调用的例子 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)); end: #p在约束Zero(a11/a12)下一定为０,ord是参变元 IsZeroUnderCons:=proc(a11,a12,p,ord) option remember; RETURN(IsEmpty(a11,{op(a12),p},ord)); end: #约束多项式conspolyset形式:[[{a11}，{a12}]，[{a21}，{a22}]] UnAmbiguous:=proc(conspolyset,var_para::list,term_order) local p,vp,result,a11,a12,a21,a22,newa11,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}; if nops(a21) <>0 then p:=numer(normalf(p,a21,term_order(op(vars)))); fi; AAA if (_field_characteristic_ <>0 ) then p :=expand(p) mod _field_characteristic_ fi; if nops(a11) <>0 then p:=numer(normalf(p,a11,term_order(op(paras)))); fi; if (_field_characteristic_ <>0 ) then p :=expand(p) mod _field_characteristic_ 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; #lprint("p="); #lprint(p); lc:=leadcoeff(p,term_order(op(vars))); #lprint("lc="); #lprint(lc); 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); newa11:=CutPolyByFactor(a11,lc); result:=UnAmbiguous([[{op(newa11),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); newa11:=CutPolyByFactor(a11,lc); result:=UnAmbiguous([[{op(newa11),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 do if divide(rt,i,'q') then rt:=q; fi; od; rt; end: <1>list1中的多项式已不可约 <2>去掉pol中在list1里列出的因子且忽略重数 CutNozeroFactors:=proc(pol,list1) local i,q,sqf_pol,rt; sqf_pol:=sqrfree(pol); sqf_pol:=sqf_pol[2]; for i from 1 to nops(sqf_pol) do sqf_pol[i][1]:=DivideList(sqf_pol[i][1],list1); od; rt:=1; for i from 1 to nops(sqf_pol) do rt:=rt*(sqf_pol[i][1]); od; rt; end: 去掉list1中以pol为因子的多项式 CutPolyByFactor:=proc(list1,pol) local p,rt,q; rt:={}; for p in list1 do if not divide(p,pol,'q') then rt:={op(rt),p};fi; od; rt; end: <1>通过gb对等于零的约束部分做normalform来化简 <2>因为gb的首系数一定不为零，通过对gb做primpart来化简 SimplifyConsgb:=proc(consgb,var_para,term_order) local rt,i,temp,paras,vars; vars:=var_para[1]; paras:=var_para[2]; rt:={}; for i in consgb do temp:=map(normalf,i[2],i[1][1],term_order(op(paras))); temp:=map(primpart,map(numer,temp),vars); rt:={op(rt),[i[1],temp]}; od; rt; end: FactorList:=proc(pol) local i,list1,list2; list1 := factors(pol)[2]; list2 := {}; for i in list1 do list2 := {op(list2),i[1]}; od; list2 end: 如果p在constrain的条件下为0在返回1,np=0 如果p在constrain的条件下为非0的常数返回10,np=10 如果p在constrain的条件下首项系数一定不是0则返回-1 其他情形返回0。np=pol #paracons=[ [part=0 ],[part <>0]] 是参数的约束; nonzerolc:=proc(pol,paracons,var_para,np,term_order) local lc,vars,paras,cons1,cons2,np0,t,pol2; global _field_characteristic_ ; vars:=var_para[1]; paras:=var_para[2]; cons1:=paracons[1]; cons2:=paracons[2]; lc:=leadcoeff(pol,term_order(op(vars))); if pol=0 then np:=0; RETURN(1); fi; if type(pol,'constant') then np:=10; RETURN(10); fi; if type(lc,'constant') then np:=pol; RETURN(-1); fi; lc在约束下一定为０ if IsZeroUnderCons(cons1,cons2,lc,paras) then pol2:=reductum(pol,term_order(op(vars))); if (_field_characteristic_ <>0 ) then pol2:=pol2 mod _field_characteristic_ fi; t:=nonzerolc(pol2,paracons,var_para,'np0',term_order); np:=np0; RETURN(t); lc在约束下一定不为０ elif NotZeroUnderCons(cons1,cons2,lc,paras) then if lc=pol then np:=10; RETURN(10); else np:=pol; RETURN(-1); fi; else np:=pol; RETURN(0); fi; end: reductum:=proc(p,ord) local res; expand(p-leadcoeff(p,ord)*leadterm(p,ord)); end: gb:=proc(ps,var_para,term_order) local res; res:=gb0([[{},{}],[{},{op(ps)}]],var_para,term_order); res:=SimplifyConsgb(res,var_para,term_order); end: gb0:=proc(conspolyset,var_para,term_order) local i,branches, res,temp; res:={}; branches:=UnAmbiguous(conspolyset,var_para,term_order); for i in branches do temp:=gb1(i,var_para,term_order); res:=res union temp; od; res; end: Step 1 gb1:=proc(conspolyset,var_para,term_order) local i,j,T,paracons,pset,pset0,pset1,pol,pset3,l,nonzerosp,sp,bs2,vars,newpol,cr,nf,np; global _field_characteristic_; T:={}; vars:=var_para[1]; paracons:=conspolyset[1]; pset:=conspolyset[2][1]; pset0:=pset; pset1:={}; if nops(pset)=0 then RETURN({[paracons,[0]]}) fi; if nops(pset) =1 then RETURN({[paracons,[op(pset)]]}) fi; Step 2----- Normalization------ for i in pset do pset0:=pset0 minus {i}; nf:=normalf(i,pset0 union pset1,term_order(op(vars))); pol:=expand(numer(nf )); cr:=nonzerolc(pol,paracons,var_para,'newpol',term_order); newpol:=numer(normalf(newpol,pset0 union pset1,term_order(op(vars)))); if (_field_characteristic_ <>0 ) then newpol :=expand(newpol) mod _field_characteristic_ fi; #lprint("newpol=");#lprint(newpol); if cr=1 then ; elif cr=10 then RETURN( { [paracons,[1]] }); elif cr=-1 then if newpol <>0 then pset1:=pset1 union {newpol} fi; else pset3:=[paracons,[ pset0 union pset1 ,{newpol} ]]; RETURN( gb0(pset3,var_para,term_order) ); fi; od; Step 3----- S-polynomial------# l:=nops(pset1); if l=1 then RETURN( {[paracons,[op(pset1)]] }) fi; nonzerosp:={}; for i to l do for j from i+1 to l do sp:=expand(numer(normalf(spoly(pset1[i],pset1[j],term_order(op(vars))),pset1,term_order(op(vars))))); if( _field_characteristic_ <>0) then sp:= expand(sp) mod _field_characteristic_ fi; #lprint("sp="); #lprint(sp); if (nonzerolc(sp,paracons,var_para,'np',term_order)<>1) then nonzerosp:=nonzerosp union {np} fi; od; od; if nonzerosp={} then RETURN( {[paracons,inter_reduce([op(pset1)],term_order(op(vars)) )] } ) ; else bs2:=[paracons,[{op(pset1)},nonzerosp ]]; T:=T union gb0(bs2,var_para,term_order); fi; T: end: # solution number # basis:=proc(a,b,ord) local i,j,k1,k;k:={}; for i in a do k1:=normalf(i,b,ord); if(k1<>0) then k:=k union {k1}; end if; end do; return k; end: mulset:=proc(a,b) local i,j,k; k:={}; if b={} then return a; fi; for i in a do for j in b do k:=k union {i*j}; od; od; return k; end: create:=proc(a,b,num) local i,j,k,k1,b1; b1:=b union {1}; k1:=b1 union mulset(a,b1); if num=0 then return k1;fi; k:=create(a,k1,num-1); end: havefinitesolution:=proc(l,var) local i,lvs,lv1; lv1:={}; for i in l do lvs:=indets(i) intersect {op(var)}; if nops(lvs)=1 then lv1:=lv1 union lvs; fi; od; if ({op(var)} minus lv1 )={} then return true; else return false; fi; end: comp:=proc(a,b) local i,j,k,l;k:={}; for i in b do for j in a do if j=i then k:=k union {j}; fi; od; od; return k; end: getmostdegree:=proc(l) local i,j,k1,k2; k2:={}; k1:={}; for i in l do k1:={op(i)}; for j in k1 do k2:=k2 union {op(j)}; od; od; return k2; end: ltt:=proc(l,vorder) local i,k; k:={}; for i to nops(l) do k:={leadterm(l[i],vorder)} union k; end do; return k; end: GetDim:=proc(lts,vars) local i,branch,l,rt; rt:=nops(lts); branch:=Nrs(lts,[op(vars)],{}); branch:=[op(branch)]; for i to nops(branch) do l:=nops(branch[i]); if l<rt then rt:=l;fi; od; return nops(vars)-rt; end: for reduced groebner basis sumofdegree:=proc(lts) local i,rt; rt:=0; for i in lts do if nops(indets(i))=1 then rt:=rt+degree(i); fi; od; return rt; end: solution1:=proc(ps,ord,term_order) local i,j,k,cgb,lt,mosttotaldeg,dim; cgb:=gb(ps,ord,term_order); cgb:=[op(cgb)]; for i to nops(cgb) do lt:=ltt(cgb[i][2],term_order(op(ord[1]))); if havefinitesolution(lt,ord[1]) then mosttotaldeg:=sumofdegree(lt); k:=basis(create({op(ord[1])},{},mosttotaldeg),lt,term_order(op(ord[1]))); cgb[i]:=[op(cgb[i]),k,cat("finite solution with number: ", nops(k))]; elif cgb[i][2]=[1] then cgb[i]:=[op(cgb[i]),"no solution"]; elif cgb[i][2]=[0] then cgb[i]:=[op(cgb[i]),cat("infinite solution with dimension: ",nops(ord[1]))]; else dim:=GetDim(lt,ord[1]); cgb[i]:=[op(cgb[i]),[lt],cat("infinite solution with dimension: ", dim)]; fi; od; return {op(cgb)}; end: 第4个参数可省,如果有,必须加单引号! solution:=proc(ps,ord,term_order,pgb) local i,j,k,cgb,lt,mosttotaldeg,dim; cgb:=gb(ps,ord,term_order); cgb:=[op(cgb)]; if nargs=4 then pgb:={op(cgb)};fi; for i to nops(cgb) do lt:=ltt(cgb[i][2],term_order(op(ord[1]))); if havefinitesolution(lt,ord[1]) then mosttotaldeg:=sumofdegree(lt); k:=basis(create({op(ord[1])},{},mosttotaldeg),lt,term_order(op(ord[1]))); cgb[i][2]:=cat("finite solution with number: ", nops(k)); elif cgb[i][2]=[1] then cgb[i][2]:="no solution"; elif cgb[i][2]=[0] then cgb[i][2]:=cat("infinite solution with dimension: ",nops(ord[1])); else dim:=GetDim(lt,ord[1]); cgb[i][2]:= cat("infinite solution with dimension: ", dim); fi; od; {op(cgb)}; end: FindRightGB:=proc(paragb) local gb,pgb,rt,i,j; j:=0; pgb:=[op(paragb)]; for i to nops(pgb) do gb:=pgb[i]; if nops(gb[1][1])=0 then gb[1][1]:={0} fi; if nops(gb[1][2])=0 then gb[1][2]:={1} fi; if {op(gb[1][1])}={0} and (gb[1][2] intersect {0})<>{0} then j:=j+1; rt:=gb[2]; fi; od; if j<>1 then print("ambiguous occur!");fi; rt; end: EvalAll:=proc(op,vars,vals) local i,rt; rt:=op; for i to nops(vars) do rt:=eval(rt,vars[i]=vals[i]); od; rt; end: #把ps中的多项式设成首项系数为１ SetLC1:=proc(ps) local i,p,rt; rt:={}; for i to nops(ps) do p:=expand(ps[i]); if p=0 then rt:={op(rt),p}; else p:=p/lcoeff(p); rt:={op(rt),p}; fi; od; rt; end: 检验约束groebnerbasis程序的正确性 check:=proc(paravalue,paragb,ps,ord,term_order) local i,gb,gb1,spec_paragb,spec_ps; if nops(paravalue)<>nops(ord[2]) then return "Error in parameter assignment: parameter number is NOT matched!" fi; spec_ps:=EvalAll(ps,ord[2],paravalue); spec_ps:={op(spec_ps)} minus {0}; spec_ps:=[op(spec_ps)]; gb1:=gbasis(spec_ps,term_order(op(ord[1]))); spec_paragb:=EvalAll(paragb,ord[2],paravalue); gb:=FindRightGB(spec_paragb); if SetLC1(gb)=SetLC1(gb1) then print( "two results are the same polynomial set!"); else print("two results are NOT the same polynomial set!"); fi; print("greobner basis gotten by parameter assignment at begin:"); print(gb1); print("greobner basis gotten by CGB method:"); gb; end: # 例子 # #参数： #ps:=[a*x^2-b*x,b*y^2-c*x,c*z^2-a*x+b]; #ord:=[[x,y,z],[a,b,c]]; #termord:=tdeg; #'pgb';　是一个变量名，该参数可缺省 #调用１： #gb(ps,ord,termord) #返回： #ps的约束GroebnerBasis #调用２： #solution(ps,ord,termord,'pgb');　 #返回： #ps的不同约束下解的个数。 #pgb为ps的约束GroebnerBasis #调用３： #solution1(ps,ord,termord);　 #返回： #ps的约束GroebnerBasis及相应解的个数（若有限，还列出不属于首理想的单项式）。 #参数： #paraval:=[1,2,3];　是参数[a,b,c]的特定值 #pgb; ps的约束GroebnerBasis #调用4： # check(paraval,pgb,ps,ord,termord); #该函数比较pgb在特定参数值paraval下的GroebnerBasis与先对参数取值后再求的GroebnerBasis是否相同。 #同时也列出这两个GroebnerBasis。此函数可做为全文的测试函数。 #返回： #pgb在特定参数值paraval下的GroebnerBasis computing partioned-parametric Grobner basis 2005.11.4 END END END END END END END END END END END END END END for computing the projection of quasi-variety 2005.11.4# read "wsolve.txt" firstly if not type(wsolve,procedure) then #lprint(" Warning: please read 'wsolve.txt' firstly ! "): fi: -----main function------------------------------------------------ Project(ps,ds,ord,pord) //old main function Sproject(result_of_project,pord) //no component is redundant Example: [[ASC,[a,b]],...] == (ASC=0 and a*b<>0 ) or ... Proj(ps,ds,ord,pord) //new main function Sproj(result_of_project,pord) //no component is redundant Example: [[ASC,[a,b]],...] == (ASC=0 and ( a<>0 or b<>0) ) or ... ------------------------------------------------------------------- #-----------Part1------<algebraic case>--------------- # project(ps,ds,ord,pord) # ord:= all variables (descending order list) # pord:=variables to eliminate (as above) #----------------------------------------------------- #QV (QuasiVariety) ; ASC (ascendant characteristic set) # get the product of the initials of cs 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: GetIP:=proc(cs,ord) local i,rt; rt:=1; for i from 1 to nops(cs) do rt:=rt*Initial(cs[i],ord); od; rt; end: # get the product of the polys in ds. GetProduct:=proc(ds) local jds,d; jds:=1; for d in ds do jds:=jds*d; od; jds; end: #get the union of two lists, denoted as list also. ListUnion:=proc(list1,list2) local rt; rt:={op(list1)} union {op(list2)}; rt:=[op(rt)]; end: # inverse a list Inverse:=proc(lis) local i,result; result:=[]; for i from nops(lis) to 1 by -1 do result:=[op(result),lis[i]]; od; result; end: # inverse every list in a list. InverseAll:=proc(lislis) local rt,l; rt:=[]; for l in lislis do rt:=[op(rt),Inverse(l)]; od; end: # eliminate such ASC in cslist that Premas(jds,ASC,ord)=0 Newcslist:=proc(cslist,jds,ord) local rt,cs; rt:=[]; for cs in cslist do if Premas(jds,cs,ord)<>0 then rt:=[op(rt),cs] fi; od; rt; end: # factors of a polynomial w.r.t ord; # facts(x^2*y+x*y-x,[x,y]); [x,x*y+y-1] # facts(x^2*y+x*y-x,[y]); [x*y+y-1] # facts(x^2*y+x*y-x,[u,v,w]); [] Ordfacts:=proc(pol,ord) local i,temp,facl,result; if type(pol,polynom) then result:=[]; facl:=factors(pol); facl:=facl[2]; if facl=[] then RETURN([]) fi; for i from 1 to nops(facl) do temp:=facl[i]; if dClass(temp[1],ord)<>0 then result:=[op(result),temp[1]]; fi; od; RETURN(result); else #lprint(`Warning: facts is called with no polynomial`); RETURN(pol); fi; end: # put p to proper position in notzerolist s.t. MV(ASC[i-1])<=MV(p)<MV(ASC[i]) PutToNotZero:=proc(p,notzerolist,mvl,ord) local i,cp,rt; rt:=notzerolist; cp:=Class(p,ord); if cp=0 then RETURN(rt);fi; for i from 1 to nops(mvl) do if cp<mvl[i] then rt:=subsop(i=rt[i]*p,rt);break; fi; od; i; if i=nops(mvl)+1 then rt:=subsop(-1=rt[-1]*p,rt) ; fi; rt; end: # all factors(noted fpolys) of polys in ds and initials, # notzerolist[i]= the product of p, # s.t. <1> p in fpolys and <2> MV(ASC[i-1])<=MV(p)<MV(ASC[i]) CreatNotZeroList:=proc(cs,ds,mvl,ord) local i,fs,p,notzerolist; notzerolist:=[]; for i from 0 to nops(cs) do notzerolist:=[op(notzerolist),1]; od; fs:=[]; for p in ds do fs:=ListUnion(fs,Ordfacts(p,ord)); od; for p in fs do notzerolist:=PutToNotZero(p,notzerolist,mvl,ord); od; notzerolist; end: #a front part of ASC whose polynomials contain no variables in pord SubASC:=proc(asc,pord) local i, result; result:=[]; for i from 1 to nops(asc) do if Class(asc[i],pord)=0 then result:=[op(result),asc[i]]; else break; fi; od; RETURN(result); end: #judge the result of ProjectASC #if it is [[...],1] return 1; #if it is [[...],0] return 0; #else return 3; IsTrival:=proc(qv,ord) local rt; if nops(qv)<>2 then rt:=3; elif not type(qv[2],polynom) then rt:=3; elif Class(qv[2],ord)<>0 then rt:=3; elif qv[2]=0 then rt:=0; else rt:=1; fi; RETURN(rt); end: # get variables of pol which are of higher order then the highest # polynomial in asc. HigherOrdVars:=proc(pol,mvl,len,ord) local vl,rt,var,index; vl:=indets(pol) intersect {op(pord)}; rt:=[]; if len<=1 then RETURN(rt);fi; for var in vl do index:=Class(var,ord); if index>mvl[len-1] then rt:=[op(rt),var]; fi; od; rt; end: # decompose the pol in [[ASC],pol] to ['UN',[[ASC],pol1],...] # so as to eliminate variables in varl. DecompList:=proc(qv,subasc,notzerolist,varList,ord) local i,flag,lenasc,rt,pol,coeffsList,coeff; if varList=[] then RETURN(qv);fi; pol:=qv[2]; if pol=0 then rt:=[subasc,pol]; else flag:=0; lenasc:=nops(qv[1]); for i from 1 to lenasc do if notzerolist[i]!=1 then flag:=1; break;fi; od; coeffsList:=[coeffs(expand(pol),varList)]; if nops(coeffsList)>1 then rt:=["UN"]; else rt:=[]; fi; for coeff in coeffsList do if Class(coeff,ord)=0 and flag=0 then RETURN([subasc,1]); else rt:=[op(rt),[qv[1],coeff]]; fi; od; fi; if nops(rt)=1 then rt:=rt[1];fi; rt; end: #when ASC has no variables in varl . # decompose the pol in [[ASC],pol] to ['UN',[[ASC],pol1],...] # so as to eliminate variables in varl DecompList1:=proc(qv,subasc,notzerolist,varList,ord) local i,flag,lenasc,rt,pol,coeffsList,coeff; if varList=[] then RETURN(qv);fi; pol:=qv[2]; if pol=0 then rt:=[subasc,pol]; else flag:=1; lenasc:=nops(qv[1]); for i from 1 to lenasc do flag:=flag*notzerolist[i] ; od; coeffsList:=[coeffs(expand(pol),varList)]; if nops(coeffsList)>1 then rt:=["UN"]; else rt:=[]; fi; for coeff in coeffsList do if Class(coeff,ord)=0 and flag=1 then RETURN([subasc,1]); else rt:=[op(rt),[qv[1],coeff*flag]]; fi; od; fi; if nops(rt)=1 then rt:=rt[1];fi; rt; end: # eliminate the highest polynomial in ASC, within [[ASC],pol] ElimPoly:=proc(qv,notzerolist,mvl,subasc,ord,pord) local rt,newasc,asc,pol,hPol,hVar,hDeg,rem,newPol,tempList,varList,len; asc:=qv[1]; len:=nops(asc); newasc:=subsop(len=NULL,asc); pol:=qv[2]; hPol:=asc[len]; hVar:=Mainvar(hPol,ord); hDeg:=degree(hPol,hVar); if degree(pol,hVar)=0 then rt:=[newasc,pol*notzerolist[len]]; RETURN(rt); else rem:=sprem(pol^hDeg,hPol,hVar); if rem=0 then rt:=[subasc,0]; RETURN(rt); fi; newPol:=rem*notzerolist[len]; tempList:=[newasc,newPol]; varList:=HigherOrdVars(newPol,mvl,len,ord,pord); rt:=DecompList(tempList,subasc,notzerolist,varList,ord); RETURN(rt); fi; end: # project wth one QV formed by one ASC(prepair share done by ProjASC) ProjectASC:=proc(qv,notzerolist,mvl,subasc,lensub,ord,pord) local pol,asc,nproj,tempList,sign,rt,nozero,i; pol:=qv[2]; asc:=qv[1]; # get the nonzero polynomial to prject with here. nozero:=pol; for i from 1 to nops(asc) do nozero:=nozero*notzerolist[i]; od; # recursion ending conditions if nops(asc)=0 or nops(asc)=lensub then if Class(nozero,pord)>0 then rt:=DecompList1(qv,subasc,notzerolist,pord,ord); elif Class(nozero,ord)>0 then rt:=[asc,nozero]; else rt:=[asc,Get1or0(nozero)]; fi; RETURN(rt); elif Class(nozero,pord)=0 then if Class(nozero,ord)>0 then rt:=[subasc,nozero]; else rt:=[subasc,Get1or0(nozero)]; fi; RETURN(rt); fi; # eliminate the highest polynomial in asc of QV. tempList:=ElimPoly(qv,notzerolist,mvl,subasc,ord,pord); # whether have many branchs in the result of ElimPoly. if not type(tempList[1],string) then rt:=ProjectASC(tempList,notzerolist,mvl,subasc,lensub,ord,pord); else rt:=[tempList[1]]; for i from 2 to nops(tempList) do nproj:=ProjectASC(tempList[i],notzerolist,mvl,subasc,lensub,ord,pord); # deal with trival cases. sign:=IsTrival(nproj,ord); if sign=0 then next; elif sign=1 then rt:= [subasc,1];break; else rt:=[op(rt),nproj]; fi; od: # rt:=["UN"] if nops(rt)=1 then rt:=[subasc,0]; fi; # rt:=["UN",[...]] if nops(rt)=2 and type(rt[2],list) then rt:=rt[2];fi; fi; rt; end: # project with one QV formed by one ASC ProjASC:=proc(cs,ds,ord,pord) local mvl,notzerolist,p,qv,subasc,lensub,rt; # mainvars order list mvl:=[]; for p in cs do mvl:=[op(mvl),Class(p,ord)]; od; # notzerolist notzerolist:=CreatNotZeroList(cs,ds,mvl,ord); qv:=[cs,notzerolist[-1]]; subasc:=SubASC(cs,pord); lensub:=nops(subasc); rt:=ProjectASC(qv,notzerolist,mvl,subasc,lensub,ord,pord); rt:=Format1(rt); end: #whether input QuasiVariety is [[[],[1]]] QVIsTruth:=proc(qv) if qv=[[[],[1]]] then RETURN(true); else RETURN(false); fi; end: #----main func.------- # project with the QV [ps,ds], eliminate variables in pord # return [] if zero(ps) is empty # return [[[],[0]]] if zero(ps/ds) is empty # else return [ [[ASC],[p1,p2,...]] , ...] Project:=proc(ps,ds,ord,pord) local jds,cslist,cs,p,ini,qvs,newds,subasc,rt; jds:=GetProduct(ds); cslist:=wsolve(ps,ord); #printf("The characteristic sets are:"); lprint(cslist); rt:=[]; if nops(cslist)<>0 then cslist:=Newcslist(cslist,jds,ord); if nops(cslist)=0 then RETURN([]); fi; cslist:=InverseAll(cslist); for cs in cslist do #add initials to ds newds:={op(ds)}; for p in cs do ini:=Initial(p,ord); if Class(ini,ord)>0 then newds:={op(newds),ini}; fi; od; newds:=[op(newds)]; qvs:=ProjASC(cs,newds,ord,pord); subasc:=qvs[1][1]; qvs:=Simplify1(qvs,subasc,ord); if QVIsTruth(qvs) then rt:=qvs; RETURN(rt); else rt:=ListUnion(qvs,rt); fi; od; fi; if nops(rt)=0 then RETURN([[[],[0]]]); fi; rt; end: the subasc in qvs is the same one AddQVS:=proc(rt,qvs) local newrt,i,j,subasc; newrt:=rt; subasc:=qvs[1]; j:=0; for i from 1 to nops(newrt) do if subasc=newrt[i][1] then newrt[i][2]:=ListUnion(qvs[2],newrt[i][2]); j:=1; break;fi; od; if j=0 then newrt:=[op(newrt),qvs];fi; newrt; end: SqfreeList:=proc(plist) local j,p,q,rt,sq; rt:=[]; for p in plist do j:=1; sq:=sqrfree(p); sq:=sq[2]; for q in sq do j:=j*q[1]; od; rt:=[op(rt),j]; od; rt; end: SqfreeGB:=proc(plist,ord) local rt,gb1,gb2; gb1:=plist; gb2:=Groebner[gbasis](gb1,plex(op(ord))); gb2:=SqfreeList(gb2); while(gb1<>gb2) do gb1:=gb2; gb2:=Groebner[gbasis](gb1,plex(op(ord))); gb2:=SqfreeList(gb2); od; gb2; end: GBsimplify:=proc(qvs,ord) local qv,rt; rt:=[]; for qv in qvs do rt:=[op(rt),[qv[1],SqfreeGB(qv[2],ord) ] ]; od; rt; end: Congre:=proc(qvs) local rt,qv,notzero; notzero:=[]; for qv in qvs do notzero:=[op(notzero),qv[2]]; od; notzero:={op(notzero)}; notzero:=[op(notzero)]; rt:=[qvs[1][1],notzero]; end: #new project Proj:=proc(ps,ds,ord,pord) local jds,cslist,cs,p,ini,qvs,newds,subasc,rt; jds:=GetProduct(ds); cslist:=wsolve(ps,ord); #printf("The characteristic sets are:"); lprint(cslist); rt:=[]; if nops(cslist)<>0 then cslist:=Newcslist(cslist,jds,ord); if nops(cslist)=0 then RETURN("Characteristic Set is empty set"); fi; cslist:=InverseAll(cslist); for cs in cslist do #add initials to ds newds:={op(ds)}; for p in cs do ini:=Initial(p,ord); if Class(ini,ord)>0 then newds:={op(newds),ini}; fi; od; newds:=[op(newds)]; qvs:=ProjASC(cs,newds,ord,pord); subasc:=qvs[1][1]; qvs:=Simplify2(qvs,subasc,ord); if qvs=[[[],1]] then RETURN([[[],[1]]]); else qvs:=Congre(qvs); rt:=AddQVS(rt,qvs); fi; od; fi; if nops(rt)=0 then RETURN([[[],[0]]]); fi; rt:=GBsimplify(rt,ord); end: # return 1 if zero(ps,ds)=empty # else return 0 IsEmpty:=proc(ps,ds,ord) local jds,cslist,cs,p,ini,qvs,newds,subasc; if nops(ps)=0 then RETURN(false);fi; jds:=GetProduct(ds); cslist:=wsolve(ps,ord); #printf("The characteristic sets are:"); lprint(cslist); if nops(cslist)<>0 then cslist:=Newcslist(cslist,jds,ord); if nops(cslist)=0 then RETURN(true); fi; cslist:=InverseAll(cslist); for cs in cslist do #add initials to ds newds:=ds; for p in cs do ini:=Initial(p,ord); if Class(ini,ord)>0 then newds:=[op(ds),ini]; fi; od; qvs:=ProjASC(cs,newds,ord,ord); subasc:=qvs[1][1]; qvs:=Simplify1(qvs,subasc,ord); if QVIsTruth(qvs) then RETURN(false); fi; od; fi; RETURN(true) end: #return -1 if zero([ ps, [op(ds),p] ])=empty #return 1 if zero([ [op(ps),p],ds ])=empty #else return 0 whichcase:=proc(ps0,ds0,p,ord) local ps,ds; ps:=[op(ps0)]; ds:=[op(ds0)]; if IsEmpty(ps,[op(ds),p],ord,ord) then RETURN(-1) fi; if IsEmpty([op(ps),p],ds,ord,ord) then RETURN(1) fi; RETURN(0); end: #---------------- Format1(qvs)----------------- #when the input is as:["UN",...] #step1: formated all the branchs of "UN". #step2: Union all the branchs. FormatUN:=proc(qvs) local i,temp,rt; rt:=[]; for i from 2 to nops(qvs) do temp:=Format1(qvs[i]); rt:=ListUnion(rt,temp); od; rt; end: # when the input is as:[[ASC],pol] # output [[[ASC],pol]] FormatList:=proc(qv) RETURN([qv]); end: # format the result of ProjectASC() # to the form: [[[ASC],pol1],[[ASC],pol2],...]; Format1:=proc(qvs) local rt; if type(qvs[1],string) then rt:=FormatUN(qvs); else rt:=FormatList(qvs); fi; rt; end: #-----------------Simplify1(qvs)--------------- #for every qausi variety in the result of Project() call doRem1 #for empty set , return []; the result form: [[[ASC],[p1,p2,...]],...] #p1<>0 and p2<>0 and ... Simplify1:=proc(qvs,subasc,ord) local qv,rt,newqv; rt:=[]; for qv in qvs do newqv:=doRem1(qv,subasc,ord); if newqv<>[] then if member(1,newqv[2]) then RETURN([[subasc,[1]]]);fi; fi; rt:=[op(rt),newqv]; od; rt; end: # a qausi variety is denoted as: [[asc],pol] # rem=premas(pol,asc,ord) # notZero:= all factors of rem of dClass w.r.t ord zero. # return [[asc], [notZero]] doRem1:=proc(qv,subasc,ord) local pol,rem,notZero,rt; pol:=qv[2]; # modified 2003.3.21 # if nops(subasc)=0 or Class(pol,ord)=0 then # if pol=0 then rt:=[]; # else rt:=[subasc,[1]]; # fi; # RETURN(rt); # fi; rem:=Premas(pol,Inverse(subasc),ord); if rem=0 then rt:=[]; else notZero:=Factorlist_chen(rem,ord); rt:=[subasc,notZero]; fi; rt; end: #-----------------Simplify2(qvs)--------------- #for every qausi variety in the result of Project() call doRem2 #for empty set , return []; the result form: [[[ASC],[p1,p2,...]],...] #p1<>0 or p2<>0 or ... #for every qausi variety in the result of Project() call doRem #for empty set , return []; the result form: [[[ASC],[p1,p2,...]],...] Simplify2:=proc(qvs,subasc,ord) local qv,rt,newqv; rt:=[]; for qv in qvs do newqv:=doRem2(qv,subasc,ord); if newqv<>[] then if newqv[2]=1 then RETURN([[subasc,1]]);fi; fi; rt:=[op(rt),newqv]; od; rt; end: # a qausi variety is denoted as: [[asc],pol] # rem=premas(pol,asc,ord) # return [[asc], rem] doRem2:=proc(qv,subasc,ord) local pol,rem,notZero,rt; pol:=qv[2]; rem:=Premas(pol,Inverse(subasc),ord); if rem=0 then rt:=[]; else notZero:=Factorlist_chen(rem,ord); rt:=[subasc,GetProduct(notZero)]; fi; rt; end: # if 0 return 0 else return 1 Get1or0:=proc(p) local rt; if p=0 then rt:=0; else rt:=1; fi; rt; end: # p=3*(x+y)^2*(x-y)*u; ord=[x,y]; result=[x+y,x-y]; # p=2; result=[1]; p=0; result=[0] Factorlist_chen:=proc(p,ord) local i,j,fl,pair,fact,rt; if Class(p,ord)=0 then rt:=[Get1or0(p)]; else rt:=[]; fl:=factors(p); fl:=fl[2]; for pair in fl do fact:=pair[1]; if Class(fact,ord)<>0 then rt:=[op(rt),fact]; fi; od; fi; rt; end: ----------Simplify---------- ----- computation of QV---------- # judge whether or not Zero(ps/ds) is empty # judge whether or not Zero(ps1/ds1) is included in Zero(ps2/ds2) # QVInclude:=proc(ps1,ds1,ps2,ds2) #------Simplify the output of Project--------# # whether or not QV is included in the union of QVS # output: 1(yes), 0(not included) QVSInclude:=proc(QV,QVS,ord) local D,newQVS,rt,temp,i,ass,proj,ps,ds,newucs; #if QVS=[] then return 0;fi; D:=GetProduct(QV[2]); newQVS:=QVS; for i from 1 to nops(QVS) do newQVS[i]:=[op(newQVS[i][1]),GetProduct(newQVS[i][2])]; od; ass:=AllSorts(newQVS,ord); for i from 1 to nops(ass) do temp:=ass[i]; ds:=[D,op(temp[2])]; if temp[1]<>[] then ps:=[op(QV[1]),op(temp[1])]; proj:=Project(ps,ds,ord,ord); # case :not empty set if proj<>[] and proj[1][2][1]<>0 then RETURN(0);fi; else proj:=ProjASC(QV[1],ds,ord,ord); # case: not empty set if proj[1][2]<>0 then RETURN(0);fi; fi; od; RETURN(1); end: # if every QV in QVS is included in the union of the others # then expunge it. Sproject:=proc(QVS,ord) local i,rt,UQV,len; rt:=[]; len:=nops(QVS); for i from 1 to len do if i=len then UQV:=[]; fi; UQV:=QVS[i+1..len]; UQV:=[op(rt),op(UQV)]; if UQV=[] then rt:=[QVS[i]];break ;fi; if QVSInclude(QVS[i],UQV,ord)=0 then rt:=[op(rt),QVS[i]]; fi; od; rt; end: Sproj:=proc(qvs,ord) local temprt,rt, tempqvs,qv,nozerop; tempqvs:=[]; #change to old form for qv in qvs do for nozerop in qv[2] do tempqvs:=[op(tempqvs), [ qv[1],[nozerop] ] ]; od; od; temprt:=Sproject(tempqvs,ord); rt:=[]; #return to new form for qv in temprt do rt:=AddQVS(rt,qv); od; rt; end: #------Simplify the output of wsolve--------# #AllSort([[a1,a3],[c1,c2,c3]]); #output:[ [[],[c1,a1]], [[],[c2,a1]], [[c3],[a1]], [[a3],[c1]], # [[a3],[c2]], [[c3,a3],[]] ] AllSorts:=proc(lists,ord) local ls,C,ls_len,p,i,reculist,list,onesort,rt; rt:=[]; ls:=lists[1]; ls_len:=nops(lists[1]); C:=Class(ls[ls_len],ord); if nops(lists)=1 then for i from 1 to ls_len-1 do p:=[[],[ls[i]]]; rt:=[op(rt),p ]; od; if C<>0 then rt:=[op(rt),[[ls[ls_len]],[]]]; fi; RETURN(rt); fi; reculist:=AllSorts(subsop(1=NULL,lists),ord); for i from 1 to ls_len-1 do for list in reculist do onesort:=list; onesort[2]:=[op(onesort[2]),ls[i]]; rt:=[op(rt),onesort]; od; od; if C<>0 then for list in reculist do onesort:=list; onesort[1]:=[op(onesort[1]),ls[ls_len]]; rt:=[op(rt),onesort]; od; fi; rt; end: # whether or not cs is included in the union of ucs # output: 1(yes), 0(not included) CSSInclude:=proc(cs,ucs,ord) local proj,IP,ass,rt,temp,i,ps,ds,newucs; if ucs=[] then return 0;fi; IP:=GetIP(cs,ord); newucs:=ucs; for i from 1 to nops(ucs) do newucs[i]:=[op(newucs[i]),GetIP(newucs[i],ord)]; od; ass:=AllSorts(newucs,ord); for i from 1 to nops(ass) do temp:=ass[i]; ds:=[IP,op(temp[2])]; if temp[1]<>[] then ps:=[op(cs),op(temp[1])]; proj:=Project(ps,ds,ord,ord); # case :not empty set if proj<>[] and proj[1][2][1]<>0 then RETURN(0);fi; else proj:=ProjASC(cs,ds,ord,ord); # case: not empty set if proj[1][2]<>0 then RETURN(0);fi; fi; od; RETURN(1); end: # if every cs in css is included in the union of the others # then expunge it. Swsolve:=proc(css,ord) local i,rt,newcss,UCS,len; rt:=[]; newcss:=InverseAll(css); len:=nops(newcss); for i from 1 to len do UCS:=newcss[i+1..len]; UCS:=[op(rt),op(UCS)]; if CSSInclude(newcss[i],UCS,ord)=0 then rt:=[op(rt),newcss[i]]; fi; od; InverseAll(rt); end: 

然后读入该文件:

read "wsolve.txt"

就可以使用wsolve()函数了

讯享网

wsolve: A Maple Package for Solving System of Polynomial Equations

Please download wsolve if you need it.

help for wsolve :

FUNCTION: wsolve - prepare the given algebraic system for solving

CALLING SEQUENCE:
wsolve(F)
wsolve(F, X)
wsolve(F, X,nonzero)
e_val(F,X,nonzero)
charset(F,X)
od_wsolve(F,X) (zero decomposition for ordinary differential poly system)
pd_wsolve(F,X) (zero decomposition for partial differential poly system)

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(asord) 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 con-
sidered in roots; however this may not stop all such solutions.


EXAMPLES:

>currentdir(): #返回当前路径,复制wsolve.txt到当前路径

>read "wsolve.txt": #读入吴方法的代码
> F := [x^2 - 2*x*z + 5, x*y^2 + y*z^3, 3*y^2 - 8*z^3]:
> wsolve(F,[z,y,x]);


{[z, y, x^2 + 5] 

[x^2-2*x*z+5, 3*y+8*x, 3*x^6-64*x^5+45*x^4+225*x^2+375]}