share/algebra/charsets/wsolve

   1 # -*- mode: maplev; -*-\r
   2 # vim: syntax=maple\r
   3 \r
   4 \r
   5 \r
   6 ######################################################### \r
   7 #11111111111111111111111111111111111111111111111111111111 \r
   8 # \r
   9 #Part I: The basic functions \r
  10 #        Class, Leader, Initial \r
  11 ######################################################### \r
  12 \r
  13 readlib(factors): \r
  14 with(Groebner):\r
  15 _field_characteristic_:=0;\r
  16 \r
  17       \r
  18 # the class of poly f wrt variable ordering ord \r
  19 Class := proc(dp,ord) \r
  20      local i,idt; \r
  21      options remember,system; \r
  22                      idt:=dIndets(dp); \r
  23                      for i to nops(ord) do \r
  24                          if member(ord[i], idt) then RETURN(nops(ord)-i+1) fi \r
  25                      od; \r
  26                      0 \r
  27                  end: \r
  28 # the highest order of the indeterminate var in dp \r
  29 # INPUT: \r
  30 #     dp : a diff pol \r
  31 #     var: an indeterminate \r
  32 # OUTPUT: the order \r
  33 # \r
  34 dOrder:=proc(dp,var) \r
  35     local idt,i,dvar; \r
  36     options remember,system;     \r
  37     idt:=indets(dp); \r
  38     dvar:=-1;      \r
  39     for i in idt do \r
  40       if has(i,var) then \r
  41            if  dvar=-1 or type(dvar,symbol) \r
  42                or (   ( not type(i,symbol) ) and ( op(1,i) > op(1,dvar) ) ) then \r
  43                dvar:=i  ; \r
  44            fi \r
  45      fi; \r
  46     od; \r
  47     if dvar=-1 then -1 \r
  48        elif   type(dvar,symbol) then 0 \r
  49        else if nops(dvar) =1 then op(1,dvar) \r
  50                 else [op(1..nops(dvar),dvar)] \r
  51              fi; \r
  52                  fi; \r
  53     end:     \r
  54        \r
  55                          \r
  56 dIndets:=proc(dp) \r
  57     local i,indet_d,idt; \r
  58     options remember,system;      \r
  59     idt:={}; \r
  60     indet_d:=indets(dp); \r
  61      \r
  62     if POLY_MODE='Algebraic' then RETURN(indet_d) fi; \r
  63      \r
  64     for i in indet_d do \r
  65       if type(i,symbol) then idt:={op(idt),i} else idt:={op(idt),op(0,i)} fi; \r
  66     od; \r
  67       idt \r
  68     end: \r
  69 \r
  70 # the initial of poly f wrt ord \r
  71 Initial := \r
  72 \r
  73    proc(f,ord) \r
  74    options remember,system; \r
  75        if Class(f,ord) = 0 then 1 \r
  76        else lcoeff(f,Leader(f,ord)); \r
  77        fi \r
  78    end: \r
  79 \r
  80 # the separant of poly f wrt ord \r
  81 Separant := \r
  82 \r
  83    proc(f,ord) \r
  84    options remember,system; \r
  85        if Class(f,ord) = 0 then 1 \r
  86        else diff(f,Leader(f,ord)); \r
  87        fi \r
  88    end: \r
  89 \r
  90 Sept:=proc(p,ord) \r
  91 local pol;\r
  92 options remember,system;  \r
  93      pol:=Separant(p,ord); \r
  94      if Class(pol,ord)=0 then RETURN(1) fi; \r
  95      pol \r
  96 end: \r
  97 \r
  98 # the set of all nonconstant factors of separants of polys in as \r
  99     \r
 100 Septs:=proc(bset,ord) \r
 101 local i,list;\r
 102   options remember,system; \r
 103     list:={}; \r
 104     for i in bset do list:=list union Factorlist(Sept(i,ord)) od; \r
 105     for i in list do if Class(i,ord)=0 then list:=list minus {i} fi od; \r
 106     list \r
 107 end: \r
 108   \r
 109     \r
 110 \r
 111 \r
 112 Lessp:= \r
 113 \r
 114 proc(f,g,ord) \r
 115 local cf,cg,cmp,df,dg,leadf,leadg; \r
 116 options remember,system; \r
 117    cf := Class(f,ord); \r
 118    cg := Class(g,ord); \r
 119    if type(f,integer) then \r
 120          true \r
 121    elif type(g,integer) then \r
 122         false      \r
 123    elif cf = 0 then \r
 124         true     \r
 125    elif cg = 0 then \r
 126         false     \r
 127    else \r
 128         leadf:=Leader(f,ord); \r
 129         leadg:=Leader(g,ord); \r
 130         \r
 131         cmp:=dercmp(leadf,leadg,ord); \r
 132         if cmp = 1 then \r
 133             true \r
 134         elif cmp=-1 then \r
 135             false \r
 136         else         \r
 137             df := degree(f,leadf); \r
 138             dg := degree(g,leadg); \r
 139             if df < dg then \r
 140                true \r
 141             elif df = dg then \r
 142                Lessp(Initial(f,ord),Initial(g,ord),ord) \r
 143             else \r
 144                false \r
 145             fi \r
 146         fi \r
 147    fi \r
 148 end:     \r
 149 \r
 150 Least:= \r
 151 \r
 152 proc(ps,ord) \r
 153 local i,lpol; \r
 154 options remember,system; \r
 155 \r
 156   lpol:=op(1,ps); \r
 157   for i from 2 to nops( ps ) do     \r
 158     if Lessp(op(i,ps),lpol,ord) then lpol:=op(i,ps) fi; \r
 159   od; \r
 160   lpol \r
 161 end: \r
 162 \r
 163 \r
 164 #############################################4 Reduced Definitions########################3 \r
 165 \r
 166 Reducedp:= \r
 167 proc(p,q,ord,T) \r
 168 local mvar,var,idt; \r
 169 options remember,system; \r
 170 \r
 171     mvar:=Leader(q,ord); \r
 172     if type(mvar,symbol) then var:=mvar else var:=op(0,mvar) fi; \r
 173      \r
 174      if POLY_MODE='Algebraic' then \r
 175 ###############Algebraic mode#################################       \r
 176        if   nargs <4 or T='std_asc' then \r
 177       if Class(p,ord) > Class(q,ord) and degree(p,mvar) < degree(q,mvar) then 1 else 0 fi; \r
 178     elif T='wu_asc' then \r
 179       if Class(p,ord) > Class(q,ord) and degree(Initial(p,ord),mvar) < degree(q,mvar) then 1 else 0 fi;        \r
 180     elif T='gao_asc' then \r
 181     ###########################in the following case , q is a list################# \r
 182       if Class(p,ord) > Class(op(1,q),ord) and Premas(Initial(p,ord),q,ord,'std_asc') <>0 then 1 else 0 fi;      \r
 183     elif T='tri_asc' then \r
 184       if Class(p,ord) > Class(q,ord) then 1 else 0 fi;      \r
 185     fi;             \r
 186 ###############Algebraic mode#################################          \r
 187      elif POLY_MODE='Ordinary_Differential' then \r
 188 ###############ODE mode#################################           \r
 189        if nargs <4 or T='std_asc' then \r
 190         \r
 191            if Class(p,ord) > Class(q,ord) and degree(p,mvar) < degree(q,mvar) \r
 192               and dOrder(p,var) <= dOrder(q,var) then \r
 193               1 \r
 194            else \r
 195               0 \r
 196            fi; \r
 197     elif T='wu_asc' then \r
 198       if Class(p,ord) > Class(q,ord) and degree(Initial(p,ord),mvar) < degree(q,mvar) then 1 else 0 fi;        \r
 199     elif T='gao_asc' then \r
 200     ###########################in the following case , q is a list################# \r
 201       if Class(p,ord) > Class(op(1,q),ord) and Premas(Initial(p,ord),q,ord,'std_asc') <>0 then 1 else 0 fi;      \r
 202     elif T='tri_asc' then \r
 203       if Class(p,ord) > Class(q,ord) then 1 else 0 fi;                  \r
 204             \r
 205        fi; \r
 206 ###############ODE mode#################################     \r
 207      elif POLY_MODE='Partial_Differential' then         \r
 208 ###############PDE mode#################################       \r
 209        if nargs <4 or T='std_asc' then \r
 210               idt:=indets(p); \r
 211               \r
 212            if dercmp(Leader(p,ord), mvar, ord)=-1 and degree(p,mvar) < degree(q,mvar) \r
 213               and   (not member(true, map(Is_proper_derivative,idt,mvar)) )  then \r
 214               1 \r
 215            else \r
 216               0 \r
 217            fi; \r
 218     elif T='wu_asc' then \r
 219       if Class(p,ord) > Class(q,ord) and degree(Initial(p,ord),mvar) < degree(q,mvar) then 1 else 0 fi;        \r
 220     elif T='gao_asc' then \r
 221     ###########################in the following case , q is a list################# \r
 222       if Class(p,ord) > Class(op(1,q),ord) and Premas(Initial(p,ord),q,ord,'std_asc') <>0 then 1 else 0 fi;      \r
 223     elif T='tri_asc' then \r
 224       if Class(p,ord) > Class(q,ord) then 1 else 0 fi;                  \r
 225             \r
 226        fi; \r
 227 \r
 228       \r
 229 ###############PDE mode#################################       \r
 230      fi; \r
 231 end:      \r
 232 \r
 233 ############################# \r
 234 # \r
 235 # Name:          Reducedset \r
 236 # Input:     ps:     a poly set \r
 237 #          p:     a polynomial \r
 238 # OUTPUT:     redset:={q | q is in ps, q is rReduced to p } \r
 239 ############################# \r
 240 \r
 241 Reducedset:= \r
 242 \r
 243 proc(ps,p,ord,T) \r
 244 local i, redset;   \r
 245 options remember,system; \r
 246     redset:={}; \r
 247     for i in ps do \r
 248       if Reducedp(i,p,ord,T)=1 then redset:={i,op(redset)} fi; \r
 249     od; \r
 250 redset \r
 251 end:      \r
 252 \r
 253 ############################################################################## \r
 254 ############################################################################## \r
 255 # the basic set of polyset ps \r
 256 Basicset:= \r
 257 \r
 258 proc(ps,ord,T) \r
 259 local qs,b,bs; \r
 260 options remember,system; \r
 261          if nops(ps) < 2 then [op(ps)] \r
 262          else \r
 263              qs:=ps;             \r
 264              bs:=[]; \r
 265              while (nops(qs) <>0) do \r
 266                 b:=Least(qs,ord); \r
 267                 bs:=[b,op(bs)]; \r
 268                 if T='gao_asc' then \r
 269                   qs :=Reducedset(qs,bs,ord,T);                 \r
 270                 else \r
 271                   qs :=Reducedset(qs,b,ord,T); \r
 272                 fi; \r
 273              od; \r
 274              bs                               \r
 275          fi \r
 276 end: \r
 277 \r
 278 ############################################################## \r
 279 # modified pseudo division: I1^s1...Ir^sr*uu = q*vv + r, \r
 280 #    where I1, ..., I_r are all distinct factors of lcoeff(vv,x) \r
 281 #    and s1, ..., sr are chosen to be the smallest \r
 282 ############################################################## \r
 283 NPrem := \r
 284 \r
 285    proc(uu,vv,x) \r
 286    local r,v,dr,dv,l,t,lu,lv,gcd0; \r
 287    options remember,system; \r
 288        if type(vv/x,integer) then subs(x = 0,uu) \r
 289        else \r
 290            r := expand(uu); \r
 291            dr := degree(r,x); \r
 292            v := expand(vv); \r
 293            dv := degree(v,x); \r
 294            if dv <= dr then l := coeff(v,x,dv); v := expand(v-l*x^dv) \r
 295            else l := 1 \r
 296            fi; \r
 297            while dv <= dr and r <> 0 do \r
 298 #                gcd0:=gcd(l,coeff(r,x,dr),'lu','lv'); \r
 299          gcd0:=gcd(l,coeff(r,x,dr)); \r
 300                lu:=simplify(l/gcd0); \r
 301                lv:=simplify(coeff(r,x,dr)/gcd0); \r
 302          lu:=l; \r
 303          lv:=coeff(r,x,dr); \r
 304                t := expand(x^(dr-dv)*v*lv); \r
 305                if dr = 0 then r := 0 else r := subs(x^dr = 0,r) fi; \r
 306                r := expand(lu*r)-t; \r
 307                dr := degree(r,x) \r
 308            od; \r
 309            r \r
 310        fi \r
 311    end: \r
 312     \r
 313     \r
 314 ################################################################################### \r
 315 ############################################################## \r
 316 # modified pseudo division: I1^s1...Ir^sr*uu = q*vv + r, \r
 317 #    where I1, ..., I_r are all distinct factors of lcoeff(vv,x) \r
 318 #    and s1, ..., sr are chosen to be the smallest \r
 319 ############################################################## \r
 320 NPremfinite := \r
 321 \r
 322    proc(uu,vv,x,ord) \r
 323    local r,v,dr,dv,l,t,lu,lv,gcd0; \r
 324    options remember,system; \r
 325    global _field_characteristic_;\r
 326        if type(vv/x,integer) then subs(x = 0,uu) \r
 327        else \r
 328            r := expand(uu); \r
 329            dr := degree(r,x); \r
 330            v := expand(vv); \r
 331            dv := degree(v,x); \r
 332            if dv <= dr then l := coeff(v,x,dv); v := expand(v-l*x^dv) \r
 333            else l := 1 \r
 334            fi; \r
 335            while dv <= dr and r <> 0 do \r
 336 #                gcd0:=gcd(l,coeff(r,x,dr),'lu','lv'); \r
 337          gcd0:=gcd(l,coeff(r,x,dr)); \r
 338                lu:=simplify(l/gcd0); \r
 339                lv:=simplify(coeff(r,x,dr)/gcd0); \r
 340          lu:=l; \r
 341          lv:=coeff(r,x,dr); \r
 342                t := expand(x^(dr-dv)*v*lv) mod _field_characteristic_; \r
 343                if dr = 0 then r := 0 else r := subs(x^dr = 0,r) mod _field_characteristic_ fi; \r
 344                r := expand(lu*r)-t mod _field_characteristic_; \r
 345 #               r :=finite_simplify(r,ord,_field_characteristic_);\r
 346                dr := degree(r,x) \r
 347            od; \r
 348            r \r
 349        fi \r
 350    end: \r
 351 #################################################\r
 352 finite_simplify:=proc(r,ord,p)\r
 353 local i,res;    \r
 354    res:=r;\r
 355    for i to nops(ord) do\r
 356      res:=subs(ord[i]^p=ord[i],res) mod p;\r
 357   od;\r
 358   res:\r
 359 end: \r
 360 \r
 361     \r
 362 ################################################################################### \r
 363 ################################################################################### \r
 364 # \r
 365 # pseudo remainder of poly f wrt ascending set as \r
 366 \r
 367 Premas := \r
 368 \r
 369    proc(f,as,ord,T) \r
 370    local r0,r1,asc;\r
 371    options remember,system;\r
 372    if nargs < 4 then asc:='std_asc' else asc:=T fi;      \r
 373      r0:=Premas_a(f,as,ord,asc);\r
 374    if POLY_MODE <> 'Partial_Differential' then RETURN(r0) fi;\r
 375    \r
 376      r1:=Premas_a(r0,as,ord,asc);\r
 377      while(r1 <> r0) do\r
 378         r0:=r1;\r
 379         r1:=Premas_a(r0,as,ord,asc);\r
 380      od;\r
 381      r1:\r
 382 end:     \r
 383      \r
 384 \r
 385 \r
 386 \r
 387 \r
 388 Premas_a := \r
 389 \r
 390    proc(f,as,ord,T) \r
 391    local remd,i,degenerate;\r
 392    options remember,system;  \r
 393    global  _field_characteristic_ ;  \r
 394        remd := f; \r
 395 #################ÏÂÃæ×¨ÃÅ´¦ÀíÓÐÏÞÓòµÄÇéÐÎ############\r
 396        if (_field_characteristic_ <> 0) then\r
 397          for i to nops(as)  do \r
 398 \r
 399             \r
 400            remd := NPremfinite(remd,as[i],Leader(as[i],ord),ord); \r
 401 \r
 402             \r
 403          od;        \r
 404        fi;\r
 405 \r
 406 \r
 407 #######################################################        \r
 408        if nargs <4 then \r
 409         \r
 410          for i to nops(as)  do \r
 411 \r
 412             \r
 413            remd := NewPrem(remd,as[i],Leader(as[i],ord)); \r
 414 \r
 415             \r
 416          od; \r
 417           \r
 418        elif T='std_asc' then \r
 419   \r
 420          for i to nops(as)  do \r
 421         \r
 422            remd := NewPrem(remd,as[i],Leader(as[i],ord)); \r
 423         \r
 424          od; \r
 425           \r
 426        elif T='tri_asc' then \r
 427         \r
 428          for i to nops(as) do \r
 429                if Class(remd,ord) = Class(as[i],ord) then                   \r
 430                        remd := NewPrem(remd,as[i],Leader(as[i],ord)); \r
 431                fi \r
 432          od;         \r
 433          \r
 434 \r
 435 ###########for wu ascending set#############         \r
 436        elif T='wu_asc' then \r
 437         \r
 438          for i to nops(as) do \r
 439              remd := WuPrem(remd,as[i],Leader(as[i],ord),ord,'degenerate'); \r
 440              if degenerate=1 then RETURN( Premas(remd,as,ord,T)) fi;\r
 441          od;       \r
 442 ##################################################                  \r
 443         \r
 444        elif T='gao_asc' then \r
 445         \r
 446          if Premas( Initial(remd,ord), as,ord,'std_asc' ) =0 then RETURN( Premas(remd,as,ord,'std_asc') ) fi;         \r
 447             for i to nops(as) do \r
 448                if Class(remd,ord) = Class(as[i],ord) then \r
 449                        remd := NewPrem(remd,as[i],Leader(as[i],ord));\r
 450                        if Premas( Initial(remd,ord), as,ord,'std_asc' ) =0 then RETURN( Premas(remd,as,ord,'std_asc') ) fi;  \r
 451  \r
 452                fi \r
 453             od;         \r
 454           \r
 455        fi; \r
 456         \r
 457         \r
 458        if remd <> 0 then numer(remd/lcoeff(remd)) else 0 fi \r
 459    end: \r
 460 \r
 461     \r
 462 ################################################################################### \r
 463 ################################################################################### \r
 464 \r
 465 # set of non-zero remainders of polys in ps wrt ascending set as \r
 466 Remset := \r
 467 \r
 468   proc(ps,as,ord,T) \r
 469   local i,set,pset,r;\r
 470   options remember,system; \r
 471   \r
 472      set:={}; \r
 473      pset:={op(ps)} minus {op(as)}; \r
 474       \r
 475      if POLY_MODE='Partial_Differential' then pset := pset union spolyset(as,ord) fi; \r
 476       \r
 477 #       print("pset=",pset); \r
 478       \r
 479       for i in pset do r:=Premas(i,as,ord,T); \r
 480             if r <>0  and Class(r,ord) =0 \r
 481                  then RETURN({1}) \r
 482                  else set:=set union {r} fi \r
 483                   od; \r
 484     set minus {0} \r
 485     end: \r
 486 ###############################misc procedures###########\r
 487 Reverse:=proc(list)\r
 488 local i,n,list1;\r
 489     n := nops(list); list1 := [seq(list[n-i+1],i = 1 .. n)]; list1\r
 490 end:\r
 491 \r
 492 ################ \r
 493 \r
 494 ################################################################################### \r
 495 # Factor the polynomials in Remaider set \r
 496 ################################################################################### \r
 497 \r
 498 Produce:=proc(rs,ord,test) \r
 499 global remfac;\r
 500 options remember,system;  \r
 501 local i,j,list1,list2; \r
 502    list2 := {}; \r
 503    for i in rs do \r
 504        list1 := Nonumber(Factorlist(i),ord); \r
 505        remfac := remfac union (list1 intersect test); \r
 506        list1 := list1 minus test; \r
 507        for j in list1 do \r
 508            if Class(j,ord) = 0 then list1 := list1 minus {j} fi \r
 509        od; \r
 510        list2 := list2 union {list1} \r
 511    od; \r
 512    for j in list2 do  if j = {} then RETURN({}) fi od; \r
 513    list2 \r
 514 end: \r
 515 \r
 516 Nonumber:=proc(list,ord) \r
 517 local i,ls;\r
 518   options remember,system;  \r
 519    ls := {}; \r
 520    for i in list do  if 0 < Class(i,ord) then ls := ls union {i} fi od; \r
 521    ls \r
 522 end: \r
 523 \r
 524 Factorlist:=proc(pol) \r
 525 local i,list1,list2; \r
 526      options remember,system; \r
 527    if (_field_characteristic_ <>0) then \r
 528      list1:=(Factors(pol) mod _field_characteristic_)[2]; \r
 529    else\r
 530      list1 := factors(pol)[2]; \r
 531    fi;\r
 532    list2 := {}; \r
 533    for i in list1 do  list2 := list2 union {numer(i[1])} od; \r
 534    list2 \r
 535 end: \r
 536 \r
 537 #input  two poly set set1,set2 \r
 538 #output \r
 539 Ltl:=proc(list1,list2) \r
 540 local set,i,j; \r
 541      options remember,system; \r
 542    set := {}; \r
 543    for i in list1 do  for j in list2 do  set := set union {{j,i}} od od \r
 544 end: \r
 545 \r
 546 Lltl:=proc(llist1,list2) \r
 547 local set,i,j; \r
 548      options remember,system; \r
 549    set := {}; \r
 550    for i in llist1 do \r
 551        for j in list2 do  set := set union {i union {j}} od \r
 552    od \r
 553 end: \r
 554 \r
 555 # Procedure Nrs: \r
 556 #input  a polynomial set  RS:={R1,R2,...,Rn} \r
 557 #output a set in which every element is a poly set. set:={set1,set2,...,sets} \r
 558 \r
 559 Nrs:=proc(rs,ord,test) \r
 560 local rm,rss,l1,i,j,rset; \r
 561      options remember,system; \r
 562    rss := Produce(rs,ord,test); \r
 563    if rss={} then RETURN({}) fi; \r
 564 #    print("rss=");print(rss); \r
 565    rset:=LProduct(rss); \r
 566    rset \r
 567 end: \r
 568 \r
 569 LProduct:=proc(inlist)    # ÊäÈëÊÇ¼¯ºÏµÄ¼¯ºÏ£¬Êä³öÒ²ÊÇ¼¯ºÏµÄ¼¯ºÏ \r
 570   local i,j,k,m,n,B,C,D,T; \r
 571      options remember,system;   \r
 572    B:=[]; \r
 573    for i from 1 to nops(op(1,inlist)) do \r
 574      B:=[op(B),{op(i,op(1,inlist))}]; \r
 575    end: \r
 576    for i from 2 to nops(inlist) do \r
 577       C:=[];D:=[]; \r
 578       for j from 1 to nops(B) do       \r
 579         if ((B[j] intersect op(i,inlist))<>{}) then   \r
 580            C:=[op(C),B[j]]; \r
 581          else D:=[op(D),B[j]]; \r
 582         end: \r
 583       end: \r
 584       B:=C; \r
 585       for m from 1 to nops(D) do \r
 586         T:=op(i,inlist); \r
 587         for n from 1 to nops(C) do \r
 588            if (nops(C[n] minus D[m])=1) then \r
 589             T:=T minus (C[n] minus D[m]);   \r
 590            end:   \r
 591          end: \r
 592        for k from 1 to nops(T) do \r
 593            B:= [op(B),D[m] union {op(k,T)}]; \r
 594         end : \r
 595       end:     \r
 596     end:   \r
 597     B:={op(B)}; \r
 598  end: \r
 599   \r
 600 LProduct_wdk:=proc(list1) \r
 601 local len,lls,i,j; \r
 602      options remember,system; \r
 603  len:=nops(list1); \r
 604  if len=0 then RETURN(list1) fi; \r
 605  lls:={}; \r
 606  for i in op(1,list1) do \r
 607    lls:=lls union { {i} } \r
 608  od; \r
 609  if len=1 then RETURN (lls) fi; \r
 610  for j from 2 to len do \r
 611    lls:=Lltl(lls, op(j,list1)); \r
 612  od; \r
 613  lls; \r
 614 end: \r
 615 \r
 616 Lltl:=proc(lls,ls) \r
 617 local res,i,j,l1,rm; \r
 618      options remember,system; \r
 619  res:={}; \r
 620  for i in lls do \r
 621    if i intersect ls ={} then \r
 622       for j in ls do \r
 623       res:= res union { i union {j} }; \r
 624       od; \r
 625    else \r
 626       res:=res union {i}; \r
 627    fi; \r
 628  od; \r
 629   \r
 630    l1 := nops(res); \r
 631    rm := {}; \r
 632    for i to l1-1 do \r
 633        for j from i+1 to l1 do \r
 634            if op(i,res) minus op(j,res) = {} then \r
 635                rm := rm union {op(j,res)} \r
 636            fi; \r
 637            if op(j,res) minus op(i,res) = {} then \r
 638                rm := rm union {op(i,res)} \r
 639            fi \r
 640        od \r
 641    od; \r
 642    res := res minus rm; \r
 643    res   \r
 644 end: \r
 645 \r
 646 ####################################################### \r
 647 # Aux functions: \r
 648 # \r
 649 # \r
 650 ###################################################### \r
 651 Nums:=proc(ps,ord) \r
 652 local i,n; \r
 653      options remember,system; \r
 654      n:=0; \r
 655      for i in ps do \r
 656        if Class(i,ord)=1 and POLY_MODE='Algebraic' then n:=n+1 fi ; \r
 657        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; \r
 658      od; \r
 659      n \r
 660 end: \r
 661 \r
 662 Emptyset:=proc(ds,as,ord) \r
 663 local i; \r
 664      options remember,system; \r
 665 \r
 666     if {op(ds)} intersect {op(as)} <> {} then RETURN(1) fi;\r
 667 #    for i in ds do if grobner[normalf](i,as,tdeg(op(ord)))=0 then RETURN(1) fi od; \r
 668    0; \r
 669 end: \r
 670 \r
 671 Non0Constp:=proc(ps,ord) \r
 672 local i; \r
 673      options remember,system; \r
 674    for i in ps do if Class(i,ord)=0 and i <> 0 then RETURN(1) fi od; \r
 675    0; \r
 676 end:     \r
 677     \r
 678 \r
 679 TestQ:=proc(ps,as,ord) \r
 680 global remfac; \r
 681 local i; \r
 682      options remember,system; \r
 683 \r
 684    for i in ps do  if grobner[normalf](i,as,ord) = 0 then remfac:=remfac union {i};RETURN(1) fi od; \r
 685    0 \r
 686 end: \r
 687 \r
 688 Init:=proc(p,ord) \r
 689 local pol; \r
 690      options remember,system; \r
 691      pol:=Initial(p,ord); \r
 692      if Class(pol,ord)=0 then RETURN(1) fi; \r
 693      pol \r
 694 end: \r
 695 \r
 696 JInitials:=proc(bset,ord) \r
 697 local pol, product; \r
 698      options remember,system; \r
 699      product:=1; \r
 700      for pol in bset do \r
 701     product:=product*Initial(pol,ord); \r
 702      od; \r
 703      product; \r
 704 end: \r
 705         \r
 706 \r
 707 Inits:=proc(bset,ord) \r
 708 local i,list; \r
 709      options remember,system; \r
 710     list:={}; \r
 711     for i in bset do list:=list union Factorlist(Init(i,ord)) od; \r
 712     for i in list do if Class(i,ord)=0 then list:=list minus {i} fi od; \r
 713     list \r
 714 end: \r
 715 \r
 716 #The following will be used in algebraci case ONLY!!! \r
 717 \r
 718 Inits1:=proc(bset,ord) \r
 719 local i,list; \r
 720      options remember,system; \r
 721     list:={}; \r
 722     for i in bset do if Class(Init(i,ord),ord)>1 then list:=list union Factorlist(Init(i,ord)) fi od; \r
 723 ###  remove the the factors with class <2 \r
 724     for i in list do if Class(i,ord) < 2 then list:=list minus {i} fi od; \r
 725     list \r
 726 end: \r
 727 \r
 728 ####################################################################### \r
 729 ####################################################################### \r
 730 # Compute the Characteristic set with FACTORIZATION      \r
 731 ####################################################################### \r
 732 ####################################################################### \r
 733 # \r
 734 # NAME:  Cs_a \r
 735 # INPUT: ps:     a polynomial set, Suppose each pol in ps is irreducible over Q. \r
 736 #      ord:   indeterminate ordering. if ord:=[z,y,x] means z>y>x; \r
 737 #      nzero: a polynomial set. Each pol in nzero is NOT equal to 0 \r
 738 #                       \r
 739 #      test:  a polynomial set, which is NOT equal to 0. \r
 740 #      T:     a symbol to decide to use which kind of ascending set \r
 741 #          T: r_asc, w_asc,g_asc,q_asc=t_asc \r
 742 # OUTPUT: a list of ascending set \r
 743 ##################################################################### \r
 744 \r
 745 Cs_a:=proc(ps,ord,nzero,test,T) \r
 746 global asc_set,INDEX,time_bs,time_rs, time_f,__testp__; \r
 747 local asc,i,j,rm,cset,test1,ps1,bset,rset, ltime_b,ltime_r,ltime_f; \r
 748 options remember,system; \r
 749 \r
 750    if Nums(ps,ord)>1 then RETURN({}) fi; \r
 751    rm := {}; \r
 752    cset := {}; \r
 753   if {op(ps)} intersect nzero <> {} then RETURN({}) fi; \r
 754   if POLY_MODE='Algebraic' then \r
 755        if nops(nzero) <> 0 and  Emptyset(nzero,ps,ord) = 1 then print("An empty set");RETURN({})  fi; \r
 756   \r
 757        if __testp__= 1 and nops(test)  <> 0 and TestQ(test ,ps,ord) = 1 then \r
 758            print("One factor of an initial has been found and it will be appended to the original polynomial set "); \r
 759 #           print("remfac=",remfac);\r
 760            RETURN({}) \r
 761        fi; \r
 762    fi; \r
 763   \r
 764 #      print("ps=",ps); \r
 765      ltime_b:=time(); \r
 766      bset := Basicset(ps,ord,T); \r
 767 #      print("bset=",bset); \r
 768      time_bs:=time_bs+time()-ltime_b; \r
 769      ltime_r:=time(); \r
 770                \r
 771      \r
 772      rset := Remset([op(ps)],bset,ord,T); \r
 773 #       print("rset=",rset); \r
 774      time_rs:=time_rs+time()-ltime_r; \r
 775 \r
 776   \r
 777    if rset={1} then RETURN({}) fi; \r
 778 # for i in rset do if Class(i,ord)=0 then RETURN({}) fi od; \r
 779    if rset = {} then \r
 780        INDEX:= INDEX+1; \r
 781        asc_set[INDEX] := bset; \r
 782        print(`A New Component`); \r
 783        RETURN({bset}) \r
 784    fi; \r
 785 #    for i in rset do \r
 786 #        if (Class(i,ord) = 0) and (i <> 0) then RETURN({}) fi \r
 787 #    od; \r
 788 \r
 789 ############## \r
 790     if POLY_MODE='Algebraic' and __testp__= 1 then\r
 791        test1 := test union Inits(bset,ord); \r
 792     else\r
 793        test1:=test;\r
 794     fi;\r
 795 ##############   \r
 796    \r
 797    \r
 798    ltime_f:=time(); \r
 799    rset := Nrs(rset,ord,nzero union test1); \r
 800    time_f:=time_f+time()-ltime_f; \r
 801    if rset = {} then RETURN({}) fi; \r
 802    cset:=map(Cs_a,map(`union`,rset,{op(bset)}),ord,nzero,test1,T); \r
 803    map(op,cset); \r
 804 end: \r
 805 \r
 806 ####################################################################### \r
 807 ####################################################################### \r
 808 # \r
 809 # NAME:  charset_a \r
 810 # INPUT: ps:     a polynomial set, Suppose each pol in ps is irreducible over Q. \r
 811 #      ord:   indeterminate ordering. if ord:=[z,y,x] means z>y>x; \r
 812 #      nzero: a polynomial set. Each pol in nzero is NOT equal to 0 \r
 813 #                       \r
 814 #      T:     a symbol to decide to use which kind of ascending set \r
 815 #          T: r_asc, w_asc,g_asc,q_asc=t_asc \r
 816 # OUTPUT: a list of ascending set \r
 817 ##################################################################### \r
 818 \r
 819 charset_a:=proc(ps,ord,nzero,T) \r
 820 global asc_set,INDEX,time_bs,time_rs, time_f; \r
 821 local asc,i,j,rm,cset,test1,ps1,bset,rset, ltime_b,ltime_r,ltime_f; \r
 822 options remember,system; \r
 823 #    if Nums(ps,ord)>1 then RETURN({}) fi; \r
 824    rm := {}; \r
 825    cset := {}; \r
 826   if {op(ps)} intersect nzero <> {} then RETURN({}) fi; \r
 827   \r
 828   if nops(nzero) <> 0 and Emptyset(nzero,ps,ord) = 1 then print("An empty set");RETURN({})   fi;   \r
 829 \r
 830   \r
 831 #      print("ps=",ps); \r
 832      ltime_b:=time(); \r
 833      bset := Basicset(ps,ord,T); \r
 834 #      print("bset=",bset); \r
 835      time_bs:=time_bs+time()-ltime_b; \r
 836      ltime_r:=time(); \r
 837      rset := Remset([op(ps)],bset,ord,T); \r
 838      rset := map(primpart,rset,ord);\r
 839 #      print("rset=",rset); \r
 840      time_rs:=time_rs+time()-ltime_r; \r
 841 \r
 842   \r
 843    if rset={1} or Non0Constp(rset,ord)=1 then RETURN([]) fi; \r
 844 # for i in rset do if Class(i,ord)=0 then RETURN({}) fi od; \r
 845    if rset = {} then \r
 846        INDEX:= INDEX+1; \r
 847        asc_set[INDEX] := bset; \r
 848        print(`A New Component`); \r
 849        RETURN(bset) \r
 850    fi; \r
 851 #    for i in rset do \r
 852 #        if (Class(i,ord) = 0) and (i <> 0) then RETURN({}) fi \r
 853 #    od; \r
 854 #    test1 := test union Inits(bset,ord); \r
 855    ltime_f:=time(); \r
 856 #    rset := Nrs(rset,ord,nzero union test1); \r
 857    time_f:=time_f+time()-ltime_f; \r
 858 #    if rset = {} then RETURN({}) fi; \r
 859 \r
 860    cset:=charset_a({op(bset)} union rset,ord,nzero,T); \r
 861    cset; \r
 862 end: \r
 863 \r
 864 charset_b:=proc(ps,ord,nzero,T) \r
 865 local cset,rs; \r
 866 options remember,system; \r
 867 \r
 868    cset := charset_a(ps,ord,nzero,T); \r
 869    rs:=Remset([op(ps)],cset,ord,T);\r
 870 #   rs:=map(primpart,rs,ord);\r
 871 \r
 872    if rs={} then RETURN(cset) fi; \r
 873    if rs={1} then RETURN([]) fi; \r
 874    if cset=[] then RETURN([]) fi; \r
 875 #   #lprint("kkkkkkkkkkkkk");\r
 876 #   cset:=charset_b({op(ps)} union {op(cset)} union rs,ord,nzero,T); \r
 877    while rs<> {} do\r
 878       cset:=charset_a({op(cset)} union rs, ord,nzero,T);\r
 879 #####we should consider the the following special case which  cset=[]      \r
 880       if nops(cset)=0 then RETURN(cset) fi;\r
 881       rs:=Remset(ps,cset,ord,T);\r
 882    od;\r
 883    cset \r
 884 end: \r
 885 \r
 886 CS_a:=proc(ps,ord,nzero,T) \r
 887 local  pset,cset,order,nonzero,asc; \r
 888 options remember, system; \r
 889    if nargs < 1 then ERROR(`too few arguments`) \r
 890      elif nargs>4 then ERROR(`too many arguments`) \r
 891    fi; \r
 892    if nargs<2 then order:=[op(indets(ps))] else order:=[op(ord)] fi; \r
 893    if nargs<3 then nonzero:={} else nonzero:=nzero fi; \r
 894    if nargs<4 then asc:=`` else asc:=T fi; \r
 895    cset:={}; \r
 896    pset:=Nrs(ps,order,nonzero); \r
 897    cset:=map(Cs_a,pset,order,nonzero,{},asc); \r
 898    [op(map(op,cset))]; \r
 899 end: \r
 900 \r
 901 ####################################################################### \r
 902 # \r
 903 # NAME:  Cs_b \r
 904 # INPUT: ps:     a polynomial set, Suppose each pol in ps is irreducible over Q. \r
 905 #      ord:   indeterminate ordering. if ord:=[z,y,x] means z>y>x; \r
 906 #      nzero: a polynomial set. Each pol in nzero is NOT equal to 0 \r
 907 #                       \r
 908 #      test:  a polynomial set, which is NOT equal to 0. \r
 909 #      T:     a symbol to decide to use which kind of ascending set \r
 910 #          T: r_asc, w_asc,g_asc,q_asc=t_asc \r
 911 # OUTPUT: a list of ascending set \r
 912 ##################################################################### \r
 913 \r
 914 Cs_b:=proc(ps,ord,nzero,test,T) \r
 915 local cset,cset1,i,j,rs,rs1; \r
 916 options remember,system; \r
 917 if Nums(ps,ord)>1 then RETURN({}) fi; \r
 918    rs1:={}; \r
 919    cset := Cs_a(ps,ord,nzero,test,T); \r
 920    cset1:=cset; \r
 921    for i in cset1 do rs:=Remset([op(ps)],i,ord,'std_asc'); \r
 922    \r
 923 #    for i in cset1 do rs:=Remset([op(ps)],i,ord,T); \r
 924 #                    if rs<>{} then if rs <> {1} then rs1:=rs1 union {rs union {op(i)} } fi; \r
 925 if rs<>{} then if rs <> {1} then rs:=Nrs(rs,ord,nzero union test); rs1:=rs1 union map(`union`,rs , {op(i)} ) fi; \r
 926                                    cset:=cset minus {i} fi \r
 927                   od; \r
 928    if rs1={} then RETURN(cset) fi; \r
 929    for j in rs1 do  cset:=cset union Cs_b(ps union j,ord,nzero,test,T) od; \r
 930    cset \r
 931 end: \r
 932 \r
 933 ####################################################################### \r
 934 # \r
 935 # NAME:  Cs_c \r
 936 # INPUT: ps:     a polynomial set, Suppose each pol in ps is irreducible over Q. \r
 937 #      ord:   indeterminate ordering. if ord:=[z,y,x] means z>y>x; \r
 938 #      nzero: a polynomial set. Each pol in nzero is NOT equal to 0 \r
 939 #                       \r
 940 #      test:  a polynomial set, which is NOT equal to 0. \r
 941 #      T:     a symbol to decide to use which kind of ascending set \r
 942 #          T: r_asc, w_asc,g_asc,q_asc=t_asc \r
 943 # OUTPUT: a list of ascending set \r
 944 ##################################################################### \r
 945 Cs_c:=proc(ps,ord,nzero,T) \r
 946 local cset,remf,ps1,i,nzeros; \r
 947 global remfac; \r
 948 options remember,system; \r
 949 \r
 950 if Nums(ps,ord)>1 then RETURN({}) fi; \r
 951     remfac:={}; \r
 952         cset:=Cs_b({op(ps)},ord,{op(nzero)},{},T); \r
 953     remf:=remfac minus {op(nzero)}; \r
 954 #    print(remf);\r
 955 #     print(remf); \r
 956   for i to nops(remf) do \r
 957 #     print(remf[i]); \r
 958 #     print(ps); \r
 959        ps1 := {op(ps)} union {remf[i]}; \r
 960        if i = 1 then nzeros := {op(nzero)} else nzeros := nzeros union {remf[i-1]} fi; \r
 961        cset := cset union Cs_c(ps1,ord,nzeros,T) \r
 962    od; \r
 963    cset \r
 964 end: \r
 965 \r
 966 CS:=proc(ps,ord,nzero,T) \r
 967 local i,pset,cset,order,nonzero,asc; \r
 968 options remember, system; \r
 969    if nargs < 1 then ERROR(`too few arguments`) \r
 970     elif nargs > 4 then ERROR(`too many arguments`) \r
 971    fi; \r
 972    if nargs < 2 then order:=[op(indets(ps))] else order:=[op(ord)] fi; \r
 973    if nargs < 3 then nonzero:={} else nonzero:=realfac({op(nzero)},ord) fi; \r
 974    if nargs < 4 then asc:='std_asc' else asc:=T fi; \r
 975    cset:={}; \r
 976    pset:=Nrs(ps,order,nonzero); \r
 977    for i to nops(pset) do cset:=cset union Cs_c(pset[i],order,nonzero,asc) od; \r
 978    [op(cset)] \r
 979 end: \r
 980 \r
 981 checknum:=0: \r
 982 Css:=proc(ps,ord,nzero,T) \r
 983 global remfac,checknum; \r
 984 local cset1,cset,i,j,nzero1, ps1,ps2,Is,Ninits; \r
 985 options remember,system; \r
 986     checknum:=checknum+1; \r
 987    remfac := {}; \r
 988    cset1 := Cs_c(ps,ord,nzero,T); \r
 989    cset := cset1; \r
 990 #          nn:=0;\r
 991       for i in cset1 do \r
 992 #         print(nops(cset1));\r
 993 #         nn:=nn+1;\r
 994 #         print(nn);\r
 995        if Class(i[nops(i)],ord)=1 and POLY_MODE='Algebraic' then Is:=Inits1(i,ord) else Is := Inits(i,ord) fi; \r
 996        if POLY_MODE='Partial_Differential' or POLY_MODE='Ordinary_Differential' then Is := Is union Septs(i,ord) fi; \r
 997        Ninits := Is minus {op(nzero)}; \r
 998 #        print("checknum=",checknum, "Ninits=",Ninits); \r
 999        for j to nops(Ninits) do \r
1000            ps1 := ({op(ps)} union {op(i)}) union {Ninits[j]}; \r
1001            if j = 1 then nzero1 := {op(nzero)} \r
1002            else nzero1 := nzero1 union {Ninits[j-1]} \r
1003            fi; \r
1004            cset := cset union Css(ps1,ord,nzero1,T) \r
1005        od \r
1006    od; \r
1007    cset \r
1008 end: \r
1009 \r
1010 realfac:=proc(ps,ord) \r
1011 local  s,res,i; \r
1012      options remember,system; \r
1013       s:=Produce(ps,ord,{}); \r
1014       res:={}; \r
1015       for i in s do \r
1016          res:=res union i; \r
1017       od; \r
1018       res; \r
1019 end:               \r
1020   \r
1021 Degenerate:=proc(ds,as,ord) \r
1022 local i,r; \r
1023      options remember,system; \r
1024        for i in ds do \r
1025          r:=Premas(i,as,ord,'std_asc'); \r
1026          if r =0 then return 1 fi; \r
1027        od; \r
1028          0; \r
1029 end: \r
1030           \r
1031 ###  POLY_MODE:= one of ['Algebraic','Ordinary_Differential','Partial_Differential']\r
1032 ### depend_relation  is like : [[x,y],[t]];          \r
1033 #### T:= one of ['std_asc','norm_asc','wu_asc','gao_asc','tri_asc']          \r
1034 CSS:=proc(ps,ord,nzero,T) \r
1035 global factor_time,prem_time, t_time,time_bs,time_rs,time_f; \r
1036 local asc,i,pset,cset,nonzero,order,result,wm; \r
1037 options remember,system; \r
1038    if nargs < 2 then order:=[op(indets(ps))] else order:=[op(ord)] fi; \r
1039    if nargs < 3 then nonzero:={} else nonzero:=realfac({op(nzero)},ord) fi; \r
1040    if nargs < 4 then asc:='std_asc' else asc:=T fi; \r
1041    if asc='norm_asc' then return Dec({op(ps)},order,nonzero,'std_asc');fi; \r
1042    cset:={}; \r
1043    factor_time:=0; \r
1044    prem_time:=0; \r
1045    time_bs:=0; \r
1046    time_rs:=0; \r
1047    time_f:=0; \r
1048     \r
1049    t_time:=time(); \r
1050    pset:=Nrs(ps,order,nonzero); \r
1051 \r
1052    \r
1053    if asc='norm_asc' then \r
1054         for i to nops(pset) do cset:=cset union Dec(pset[i],order,nonzero,'std_asc') od;\r
1055    else         \r
1056         for i to nops(pset) do cset:=cset union Css(pset[i],order,nonzero,asc) od; \r
1057    fi;\r
1058    result:=[]; \r
1059    for i in cset do if Degenerate(nonzero,i,order)<>1 then result:=[op(result),i]; fi od; \r
1060        #lprint("Timing","Total", time()-t_time, "Factoring", factor_time,"Prem",prem_time); \r
1061        #lprint("BasicSet",time_bs,"RS",time_rs,"factor",time_f); \r
1062    result; \r
1063 \r
1064 end: \r
1065 \r
1066 charset:=proc(ps,ord,nzero,T) \r
1067 global factor_time,prem_time, t_time,time_bs,time_rs,time_f; \r
1068 local asc,nonzero,order,result,cm; \r
1069 options remember,system; \r
1070    if nargs < 2 then order:=[op(indets(ps))] else order:=[op(ord)] fi; \r
1071    if nargs < 3 then nonzero:={} else nonzero:=realfac({op(nzero)},ord) fi; \r
1072    if nargs < 4 then asc:='std_asc' else asc:=T fi; \r
1073 \r
1074    prem_time:=0; \r
1075    time_bs:=0; \r
1076     \r
1077   \r
1078         \r
1079    t_time:=time(); \r
1080    \r
1081   result:=charset_b({op(ps)}, order,nonzero,asc); \r
1082 #   result:=charset_b(map(primpart,{op(ps)},order),order,nonzero,asc); \r
1083 \r
1084        #lprint("Timing","Total", time()-t_time, "BasicSet",time_bs,"Prem",prem_time); \r
1085 \r
1086    result; \r
1087 \r
1088 end: \r
1089 \r
1090 Css1:=proc(ps,ord,test,type) \r
1091 local  cset,i,j,septs,cset1,nonzero,vv; \r
1092      options remember,system; \r
1093        cset:=CSS(ps,ord,test,type); \r
1094 #       septs:={}; \r
1095        cset1:=cset; \r
1096        nonzero:={op(test)};\r
1097        for i in cset1 do \r
1098                septs:=Septs(i,ord) minus Inits(i,ord); \r
1099                if septs<>{} then\r
1100 \r
1101                        for j in septs do \r
1102                        vv:=  {op(Css1(ps union {op(i)} union {j},ord,nonzero,type))};   \r
1103                        nonzero:={op(nonzero),j};                   \r
1104                        cset:={op(cset)} union vv ;\r
1105                        od \r
1106                fi \r
1107        od; \r
1108    [op(cset)] \r
1109 end: \r
1110 \r
1111 Gsolve:=proc(ps,ord,test) \r
1112     gsolve(ps,ord,test) \r
1113 end: \r
1114 \r
1115 e_val:=proc(ps,ord,test,T) \r
1116 global factor_time,prem_time, t_time,time_bs,time_rs,time_f; \r
1117 local asc,nonzero,order,result,cm; \r
1118 options remember,system; \r
1119    if nargs < 2 then error("The input should have at least 2 parameters") fi; \r
1120    if nargs < 3 then nonzero:={} else nonzero:=realfac({op(nzero)},ord) fi; \r
1121    if nargs < 4 then asc:='std_asc' else asc:=T fi; \r
1122 \r
1123    prem_time:=0; \r
1124    time_bs:=0;               \r
1125    t_time:=time(); \r
1126    nonzero:=test;\r
1127   result:=Css1({op(ps)}, ord,nonzero,asc); \r
1128        #lprint("Timing","Total", time()-t_time, "BasicSet",time_bs,"Prem",prem_time); \r
1129 \r
1130    result; \r
1131 end: \r
1132 \r
1133  \r
1134 \r
1135 wsolve:=proc() \r
1136 local rt; \r
1137 global POLY_MODE; \r
1138 POLY_MODE:='Algebraic'; \r
1139  rt:=CSS(args); \r
1140  rt; \r
1141 end: \r
1142 \r
1143 od_wsolve:=proc() \r
1144 local rt; \r
1145 global POLY_MODE,depend_relation; \r
1146 POLY_MODE:='Ordinary_Differential'; \r
1147  if type(depend_relation,'list') and nops(depend_relation) >1 and \r
1148     type(depend_relation[2],'list') and nops(depend_relation[2])=1 then \r
1149      rt:=CSS(args); \r
1150   else \r
1151      #lprint("Please set depend_relation:=[[x,y,...],[t]] first, if x,y,... depend on t"); \r
1152      rt:=0; \r
1153   fi; \r
1154  rt; \r
1155 end: \r
1156 \r
1157 pd_wsolve:=proc() \r
1158 local rt; \r
1159 global POLY_MODE,depend_relation; \r
1160 POLY_MODE:='Partial_Differential'; \r
1161  if type(depend_relation,'list') and nops(depend_relation) >1 and \r
1162     type(depend_relation[2],'list') and nops(depend_relation[2])>1 then \r
1163      rt:=CSS(args); \r
1164   else \r
1165      #lprint("Please set depend_relation:=[[x1,x2,...],[t1,t2,...]] first, if x1,x2,... depend on t1,t2,..."); \r
1166      rt:=0; \r
1167   fi; \r
1168  rt; \r
1169 end: \r
1170 \r
1171 INDEX:=0: \r
1172 __testp__:=0:\r
1173 remfac:={}: \r
1174 prem_time:=0: \r
1175 factor_time:=0: \r
1176 time_bs:=0: \r
1177 time_rs:=0: \r
1178 time_f:=0: \r
1179 POLY_MODE:=Algebraic: \r
1180 RANKING_MODE:='cls_ord_var': \r
1181 \r
1182 # Test Program for Differentiation \r
1183 \r
1184 \r
1185 #---------------------------------------- \r
1186 # Global dependence \r
1187 #---------------------------------------- \r
1188 \r
1189 ######################################################### \r
1190 #11111111111111111111111111111111111111111111111111111111 \r
1191 # \r
1192 #Part IV: The basic functions for differential poly \r
1193 #         \r
1194 ######################################################### \r
1195 \r
1196 #Diff_Var:=[u,v,w]; # Declare variables; \r
1197 \r
1198 #Diff_Indets:=[X1,X2,X3]; # Declare Diff ndeterminates; \r
1199 \r
1200 #Lvar:=nops(Diff_Var); \r
1201 \r
1202 #---------------------------------------------------------- \r
1203 # main function \r
1204 #---------------------------------------------------------- \r
1205 \r
1206 dDiff:=proc(dp,var,n) \r
1207    local df,i; \r
1208      options remember,system;    \r
1209    df:=dp; \r
1210    if nargs=2 then RETURN (DF(dp,var)) fi; \r
1211    for i to n do \r
1212       df:=DF(df,var); \r
1213    od; \r
1214    df; \r
1215 end: \r
1216           \r
1217 #========================================================== \r
1218 #         \r
1219 #          dq  <- DF(dp, var) \r
1220 # \r
1221 #          (differetiating a diff polynomial) \r
1222 # \r
1223 # \r
1224 #          Input: dp, a diff polynomial; \r
1225 #                 var,  the variable w.r.t which to be differentiate   \r
1226 # \r
1227 #          Output: dq,  the result after differentiating dp w.r.t var; \r
1228 # \r
1229 #=========================================================== \r
1230 \r
1231 DF:=proc(dp, var) \r
1232    local dq, dp1, dp2; \r
1233      options remember,system;    \r
1234 \r
1235 #Step1 [Recursive base] \r
1236 \r
1237       if op(0,dp) <> `+` then \r
1238          dq:=DF_prod(dp,var); \r
1239          RETURN(dq); \r
1240       fi; \r
1241     \r
1242 #Step2 [Recursion] \r
1243 \r
1244      dp1:=op(1,dp); \r
1245      dp2:=normal(dp-dp1); \r
1246      dq:=normal(DF_prod(dp1,var)+DF(dp2,var));             \r
1247           \r
1248 #Step3 [Prepare for return] \r
1249 \r
1250      RETURN(dq); \r
1251       \r
1252 end: \r
1253 \r
1254 #------------------------------------------ \r
1255 # subroutines \r
1256 #------------------------------------------ \r
1257 \r
1258 #========================================================= \r
1259 # \r
1260 # \r
1261 #                der <- DF_indet(indet, var) \r
1262 #       \r
1263 #           (differentiate a differential indeterminante) \r
1264 # \r
1265 #   Input: indet, a differential indeterminate which is a \r
1266 #                 derivative of symobls in Diff_Indets w.r.t \r
1267 #                 some variables in Diff_Var; \r
1268 # \r
1269 #          var,   the variable w.r.t which to be differeniated \r
1270 # \r
1271 #   Output: der,  the derivative of indet w.r.t var \r
1272 # \r
1273 #=========================================================== \r
1274 \r
1275 DF_indet:=proc(indet, var) \r
1276          local der, p, i, index, j,Diff_Var,Diff_Indets,Lvar; \r
1277          global depend_relation; \r
1278      options remember,system;          \r
1279 \r
1280 #STEP0 [Initiation diffs]; \r
1281           \r
1282          Diff_Var:=depend_relation[2]; \r
1283          Diff_Indets:=depend_relation[1]; \r
1284          Lvar:=nops(Diff_Var); \r
1285           \r
1286 #STEP1 [Special cases] \r
1287           \r
1288          if  not member(var, Diff_Var, 'p') then \r
1289              der := diff(indet, var); \r
1290              RETURN(der); \r
1291          fi; \r
1292           \r
1293          if not member(indet, Diff_Indets) and \r
1294             not member(op(0,indet), Diff_Indets) then \r
1295             der := diff(indet, var); \r
1296             RETURN(der); \r
1297          fi; \r
1298           \r
1299 #STEP2 [Compute] \r
1300 \r
1301          index:=NULL; #Initialization \r
1302           \r
1303          if member(indet, Diff_Indets) then   \r
1304           \r
1305          # build a derivative from a diff indet       \r
1306           \r
1307             for i from 1 to Lvar do \r
1308                 if i = p then \r
1309                    index:=index,1; \r
1310                 else \r
1311                    index:=index,0; \r
1312                 fi; \r
1313             od; \r
1314             der:=indet[index]; \r
1315          else   \r
1316           \r
1317          # modify a derivative from a given one \r
1318           \r
1319            for i from 1 to Lvar do \r
1320                j:=op(i,indet); \r
1321                if i = p then \r
1322                   index:=index,(j+1); \r
1323                else \r
1324                   index:=index,j; \r
1325                fi; \r
1326            od; \r
1327            der:=op(0,indet)[index] \r
1328          fi; \r
1329           \r
1330 #STEP3 [Prepare for return] \r
1331           \r
1332           RETURN(der); \r
1333           \r
1334 end: \r
1335 \r
1336 #======================================================== \r
1337 # \r
1338 #           dq <- DF_power(dp, var) \r
1339 # \r
1340 #           (differentiating a power of a drivative of a diff indet) \r
1341 # \r
1342 #           input:  dp, a power of a diff indet \r
1343 #                   var, the variable w.r.t which to be differentiated \r
1344 # \r
1345 #           output: dq, the result after differentiating dp w.r.t var \r
1346 # \r
1347 #========================================================= \r
1348 \r
1349 DF_power:=proc(dp, var) \r
1350          local dq, indet, expon; \r
1351               options remember,system; \r
1352           \r
1353 #Step1 [Special cases] \r
1354 \r
1355       # constant case \r
1356       \r
1357       if dp = 1 then \r
1358          dq := 0; \r
1359          RETURN(dq); \r
1360       fi; \r
1361       \r
1362       # indeterminate case \r
1363       \r
1364       if op(0,dp) <> `^` then \r
1365          dq := DF_indet(dp, var); \r
1366          RETURN(dq); \r
1367       fi; \r
1368                 \r
1369 #Step2 [Compute] \r
1370 \r
1371       indet:=op(1,dp); \r
1372       expon:=op(2,dp); \r
1373       dq:=expon*indet^(expon-1)*DF_indet(indet,var); \r
1374     \r
1375 #Step3 [Prepare for return] \r
1376 \r
1377       RETURN(dq); \r
1378 \r
1379 end: \r
1380 \r
1381 #========================================================= \r
1382 # \r
1383 #        dq <- DF_prod(dp, var) \r
1384 # \r
1385 #        (differentiating a product of derivatives of some diff indets)     \r
1386 # \r
1387 #        input: dp, a product of derivatives of some diff indets   \r
1388 #               var,  the variable w.r.t which to be differentiated \r
1389 # \r
1390 #        output:  dq, the result after differentiating dp w.r.t var; \r
1391 # \r
1392 #========================================================== \r
1393 \r
1394 DF_prod:=proc(dp, var) \r
1395         local dq, dp1, dp2; \r
1396              options remember,system; \r
1397 \r
1398 #Step1 [Recursive base] \r
1399 \r
1400       if op(0,dp) <> `*` then \r
1401          dq := DF_power(dp, var); \r
1402          RETURN(dq); \r
1403      fi; \r
1404 \r
1405 #Step2 [Recursion] \r
1406 \r
1407      dp1:=op(1,dp); \r
1408      dp2:=normal(dp/dp1); \r
1409       \r
1410      dq:=normal(DF_power(dp1,var)*dp2+dp1*DF_prod(dp2,var)); \r
1411       \r
1412 #Step3 [Prepare for return] \r
1413 \r
1414      RETURN(dq); \r
1415       \r
1416 end: \r
1417                                         \r
1418 #=============================================================== \r
1419 # \r
1420 #            ord  <- torder(der) \r
1421 # \r
1422 #           computing the order of a derivative of a diff indet \r
1423 # \r
1424 #           input: der, a derivative of a diff indet \r
1425 # \r
1426 #           output: ord, ord is the order of der. \r
1427 # \r
1428 #================================================================= \r
1429 \r
1430 torder:=proc(der) \r
1431 local ord,i,Diff_Var,Lvar; \r
1432 global depend_relation; \r
1433      options remember,system; \r
1434 \r
1435 #STEP1 [Initialize] \r
1436 \r
1437       if type(der,symbol) then RETURN(0) fi; \r
1438 \r
1439          ord := 0; \r
1440               \r
1441 #STEP2 [Compute] \r
1442 \r
1443            \r
1444               for i from 1 to nops(der) do \r
1445                  ord := ord + op(i, der); \r
1446               od; \r
1447  ord; \r
1448 end: \r
1449                     \r
1450 \r
1451 #================================================================== \r
1452 #     \r
1453 #            lead <- Leader(dp, ranking,ord, _cls, _ord) \r
1454 #     \r
1455 #            Input: dp, a  nonzero diff poly; \r
1456 #                   ranking, a specified ranking; \r
1457 # \r
1458 #            Output: the Leader w.r.t. ranking; \r
1459 # \r
1460 #            Option: _cls, the class of the lead; \r
1461 #                   _ord,, the order of the lead; \r
1462 # \r
1463 #=================================================================== \r
1464 Mainvar:=proc(p,ord) \r
1465      options remember,system; \r
1466  Leader(p,ord); \r
1467  end: \r
1468 \r
1469 Leader := proc(dp, ord, _cls, _ord) \r
1470          local lead, cls, order, L, l,  der, clsord, i, j, cls1, ord1,Diff_Var,Lvar; \r
1471       global depend_relation; \r
1472            options remember,system;     \r
1473            \r
1474 #Step1 [Initialize] \r
1475 \r
1476          Diff_Var:=depend_relation[2]; \r
1477          Lvar:=nops(Diff_Var); \r
1478          L:=indets(dp); l:=nops(L); \r
1479          lead := 1; \r
1480          cls := 0; \r
1481          order := 0; \r
1482 \r
1483 #Step 2 [Algebraic mode] \r
1484 \r
1485         if POLY_MODE='Algebraic' then \r
1486          for i to nops(ord) do \r
1487            if member(ord[i],L) then RETURN(ord[i]) fi \r
1488          od; \r
1489         fi; \r
1490 \r
1491 #Step2 [cls_ord_var case] \r
1492 \r
1493          if RANKING_MODE = cls_ord_var then \r
1494             for i from  1 to l do \r
1495                 der := L[i]; \r
1496 \r
1497                 cls1:=Class(der,ord); \r
1498                 ord1:=torder(der); \r
1499                 \r
1500                 if cls1 > cls  or (cls1 = cls and ord1 > order) then \r
1501                    lead := der; \r
1502                    cls := cls1; \r
1503                    order := ord1 \r
1504                 else \r
1505                    if cls1 = cls and ord1 = order and ord1 > 0 then \r
1506                       for j from 1 to Lvar do \r
1507                          if op(j, lead) < op(j, der) then \r
1508                             lead := der; \r
1509                             cls := cls1; \r
1510                             order := ord1; \r
1511                             break; \r
1512                          else \r
1513                             if op(j, lead) > op(j, der) then; \r
1514                             break; \r
1515                             fi; \r
1516                          fi; \r
1517                      od; \r
1518                    fi; \r
1519                fi; \r
1520             od; \r
1521          fi; \r
1522           \r
1523 #Step3 [ord_cls_var case] \r
1524           \r
1525          if RANKING_MODE = ord_cls_var then \r
1526             for i from  1 to l do \r
1527                 der := L[i]; \r
1528                 cls1:=Class(der,ord); \r
1529                 ord1:=torder(der); \r
1530                 if ord1 > order or (ord1 = order and cls1 > cls) then \r
1531                    lead := der; \r
1532                    cls := cls1; \r
1533                    order := ord1 \r
1534                 elif ord1 = order and cls1 = cls and ord1 > 0 then \r
1535                      for j from 1 to Lvar do \r
1536                          if op(j, lead) < op(j, der) then \r
1537                             lead := der; \r
1538                             cls := cls1; \r
1539                             order := ord1; \r
1540                             break; \r
1541                          elif op(j, lead) > op(j, der) then; \r
1542                             break; \r
1543                          fi; \r
1544                      od; \r
1545                fi; \r
1546             od; \r
1547          fi; \r
1548           \r
1549 \r
1550 #STEP3 [Prepare for return] \r
1551 \r
1552           \r
1553          if nargs > 2 then \r
1554             _cls := cls; \r
1555          fi; \r
1556          if nargs > 3 then \r
1557             _ord := order; \r
1558          fi; \r
1559           \r
1560          RETURN(lead); \r
1561           \r
1562 end: \r
1563 \r
1564 \r
1565 \r
1566          \r
1567 \r
1568 \r
1569           \r
1570 depend:=proc(l1,l2) \r
1571 global depend_relation; \r
1572  depend_relation:=[l1,l2]; \r
1573  1: \r
1574  end: \r
1575 #============================================================ \r
1576 # \r
1577 #       {1,-1,0} <- dercmp(der1, der2) \r
1578 # \r
1579 #       input: der1, der2, two derivatives; \r
1580 #           \r
1581 # \r
1582 #       output: 0 if der1 = der2 \r
1583 #               1 if der1 < der2 \r
1584 #              -1 if der1 > der2 \r
1585 # \r
1586 #============================================================   \r
1587 dercmp:=proc(der1, der2,ord) \r
1588            local dp, lead; \r
1589                 options remember,system; \r
1590             \r
1591            if der1=der2 then RETURN(0) fi; \r
1592             \r
1593            dp := der1 + der2; \r
1594            lead:=Leader(dp,ord); \r
1595            if lead = der1 then \r
1596               RETURN(-1); \r
1597            else \r
1598               RETURN(1); \r
1599            fi; \r
1600             \r
1601 end: \r
1602 \r
1603 #if der1 is a derivative of der2 then return true else false\r
1604 Is_derivative:=proc(der1, der2) \r
1605               local dt1,dt2,i, l1,l2; \r
1606                    options remember,system; \r
1607               \r
1608               \r
1609                if type(der1,symbol) then dt1:=der1 else dt1:=op(0,der1) fi; \r
1610                if type(der2,symbol) then dt2:=der2 else dt2:=op(0,der2) fi; \r
1611             if dt1<>dt2 then RETURN(false) fi; \r
1612             if der1=der2 then RETURN(true) fi; \r
1613              \r
1614             if torder(der2)=0 then \r
1615                  RETURN(true) \r
1616             elif torder(der1)=0 then \r
1617                    RETURN(false) ; \r
1618            else \r
1619                 l1:=dOrder(der1,dt1); \r
1620                 l2:=dOrder(der2,dt1); \r
1621                 for i to nops(l1) do \r
1622                      if op(i,l1) < op(i,l2) then RETURN(false) fi; \r
1623                 od; \r
1624          fi; \r
1625                            \r
1626            true; \r
1627 end: \r
1628                                                   \r
1629 #if der1 is a proper derivative of der2 then return true else false\r
1630 Is_proper_derivative:=proc(der1, der2) \r
1631               local dt1,dt2,i, l1,l2; \r
1632                    options remember,system; \r
1633               \r
1634               \r
1635                if type(der1,symbol) then dt1:=der1 else dt1:=op(0,der1) fi; \r
1636                if type(der2,symbol) then dt2:=der2 else dt2:=op(0,der2) fi; \r
1637             if dt1<>dt2 then RETURN(false) fi; \r
1638             if der1=der2 then RETURN(false) fi; \r
1639              \r
1640             if torder(der2)=0 then \r
1641                  RETURN(true) \r
1642             elif torder(der1)=0 then \r
1643                    RETURN(false) ; \r
1644            else \r
1645                 l1:=dOrder(der1,dt1); \r
1646                 l2:=dOrder(der2,dt1); \r
1647                 for i to nops(l1) do \r
1648                      if op(i,l1) < op(i,l2) then RETURN(false) fi; \r
1649                 od; \r
1650          fi; \r
1651                            \r
1652            true; \r
1653 end: \r
1654 \r
1655 #return all the derivatives in dp of der\r
1656 proper_derivatives:=proc(dp,der) \r
1657     local i,idt,dets; \r
1658          options remember,system; \r
1659     idt:=indets(dp); \r
1660     dets:={}; \r
1661      \r
1662     for i in idt do \r
1663       if Is_proper_derivative(i,der) then \r
1664          dets:={i,op(dets)}; \r
1665          \r
1666       fi; \r
1667     od; \r
1668     dets;      \r
1669 end:      \r
1670 \r
1671        \r
1672      \r
1673      \r
1674 \r
1675 #sometimes, Maple's GCD is NOT valid for parameters which has indexes such a[1],a[2]; \r
1676 \r
1677 NPremO := \r
1678 \r
1679    proc(uu,vv,xx) \r
1680    local r,v,dr,dv,lcv,rtv,g,lu,lv; \r
1681     options remember,system; \r
1682      if degree(vv,xx)=0 then RETURN(0) fi; \r
1683      if type(vv/xx,integer) then coeff(uu,xx,0) \r
1684        else \r
1685            r := expand(uu); \r
1686            dr := degree(r,xx); \r
1687            v := expand(vv); \r
1688            dv := degree(v,xx); \r
1689             \r
1690              lcv:=lcoeff(v,xx); \r
1691 #            rtv:=sterm(v,xx); \r
1692            rtv:=v-expand(lcoeff(v,xx)*xx^degree(v,xx)); \r
1693 \r
1694            while dv <= dr and r <> 0 do \r
1695 \r
1696 #              g:=gcd(lcv,coeff(r,xx,dr)); \r
1697 \r
1698              lu:=lcv; \r
1699 \r
1700              lv:=coeff(r,xx,dr); \r
1701                   \r
1702                r:=expand(lu*(r-expand(lcoeff(r,xx)*xx^degree(r,xx)))-lv*rtv*xx^(dr-dv)); \r
1703                dr := degree(r,xx) \r
1704            od; \r
1705            r \r
1706        fi \r
1707    end: \r
1708   \r
1709   \r
1710 #######±äÁ¿yÓÐÁ½ÖÖÇéÐÎ.y and y[1]ËùÒÔÒªÐ¡ÐÄ. \r
1711 \r
1712   \r
1713 ODPrem:= \r
1714 \r
1715   proc(uu,vv,y) \r
1716     local u,x,dv,du,d,t,var; \r
1717     global depend_relation; \r
1718     options remember,system;         \r
1719     \r
1720        t:=depend_relation[2][1]; \r
1721        var:=depend_relation[1]; \r
1722   \r
1723 #        if dClass(uu,var)=0 then RETURN(uu) fi; \r
1724 \r
1725           u:=uu; \r
1726           if type(y,symbol) then x:=y else         x:=op(0,y)  fi;                 \r
1727           dv:=dOrder(vv,x); \r
1728           du:=dOrder(u,x); \r
1729           d:=du-dv; \r
1730           if d<0 then RETURN(uu) fi; \r
1731 #the following program is for ordinary differential case                 \r
1732           while d>0 do \r
1733 \r
1734             u:=NPremO(u,dDiff(vv,t,d), x[du]); \r
1735             du:=dOrder(u,x); \r
1736             d:=du-dv; \r
1737           od; \r
1738 #         \r
1739           if dv =0 then NPremO(u,vv,x) else NPremO(u,vv,x[dv]) fi;                           \r
1740 \r
1741             \r
1742    end:     \r
1743 \r
1744     \r
1745 NewPrem:= \r
1746 \r
1747   proc(uu,vv,y) \r
1748     local u,x,dv,du,d,t,var; \r
1749     global depend_relation;\r
1750     options remember,system;      \r
1751     \r
1752     if POLY_MODE='Algebraic' then \r
1753     \r
1754        NPrem(uu,vv,y); \r
1755         \r
1756     elif POLY_MODE='Ordinary_Differential' then         \r
1757     \r
1758     ODPrem(uu,vv,y);                      \r
1759      \r
1760     elif POLY_MODE='Partial_Differential' then \r
1761         \r
1762        PDPrem(uu,vv,y); \r
1763         \r
1764     fi; \r
1765             \r
1766    end:     \r
1767 \r
1768 \r
1769 \r
1770 Headder:=proc(idts) \r
1771 local hder,i,j; \r
1772      options remember,system; \r
1773   hder:=op(1,idts);   \r
1774   \r
1775 #   if type(hder,symbol) then var:=hder else var:=op(0,hder) fi; \r
1776         \r
1777   for i in idts do \r
1778   \r
1779     if torder(i) > torder(hder) then \r
1780         hder:=i ;           \r
1781     elif torder(i)=torder(her) then \r
1782        for j to nops(hder) do \r
1783           if op(j,i) > op(j,hder) then hder:=i;break; \r
1784           elif op(j,i) < op(j,hder) then break; \r
1785           fi; \r
1786        od; \r
1787     fi; \r
1788                    \r
1789    od; \r
1790     \r
1791   hder: \r
1792   \r
1793 end: \r
1794 \r
1795 \r
1796 PDPrem:=proc(dp, dq, y) \r
1797           local Diff_Indets,Diff_Var, var,dr, idt,l2,lder,da,l1,i,s; \r
1798           options remember,system; \r
1799 \r
1800 #Step1 [Special case] \r
1801 \r
1802           \r
1803     Diff_Indets:=depend_relation[1]; \r
1804     Diff_Var:=depend_relation[2]; \r
1805                \r
1806     if y=1 then RETURN(0) fi; \r
1807      \r
1808     if type(y,symbol) then var:=y else var:=op(0,y) fi; \r
1809      \r
1810     dr:=dp; \r
1811      \r
1812        idt:=proper_derivatives(dr,y); \r
1813        l2:=dOrder(y,var); \r
1814 \r
1815 #Step2 [Recursive base]         \r
1816 \r
1817        while  nops(idt)<>0 do \r
1818 \r
1819           \r
1820           lder:=Headder(idt); \r
1821           \r
1822           da := dq; \r
1823                     \r
1824 \r
1825           \r
1826           l1:=dOrder(lder,var); \r
1827 \r
1828           \r
1829           \r
1830           for i to nops(l1) do \r
1831           \r
1832            if l2=0 then s:=op(i,l1) else s:=op(i,l1) -op(i,l2) fi; \r
1833            \r
1834                   da := dDiff(da, Diff_Var[i],s); \r
1835             \r
1836           od; \r
1837           \r
1838           dr := NPremO(dr, da, lder); \r
1839           \r
1840        idt:=proper_derivatives(dr,y);           \r
1841           \r
1842     od; \r
1843 \r
1844          dr:=NPremO(dr,dq,y) ; \r
1845           \r
1846 #Step4 [Prepare for return] \r
1847 \r
1848     dr;           \r
1849           \r
1850 end:                 \r
1851           \r
1852   \r
1853     \r
1854 \r
1855 # \r
1856 # set:  A set; \r
1857 \r
1858 cmb:=proc(set0) \r
1859 local p,q, ls,set1; \r
1860      options remember,system; \r
1861     ls:={}; \r
1862     \r
1863     set1:=set0; \r
1864     \r
1865  for p in set0 do \r
1866     set1:=set1 minus {p}; \r
1867     for q in set1  do \r
1868         ls:={op(ls),[p,q]} \r
1869     od; \r
1870  od; \r
1871  ls; \r
1872 end: \r
1873   \r
1874 \r
1875 ######## modified Nov.6 ######## \r
1876 \r
1877           \r
1878 spolyset:=proc(as,ord) \r
1879 local res,bs,l,p,cls,i,poly; \r
1880      options remember,system; \r
1881     \r
1882    res:={}; \r
1883    bs:={op(as)}; \r
1884    l:=nops(bs); \r
1885    if l <= 1 then RETURN(res) fi; \r
1886     \r
1887    p:=op(1,bs); \r
1888    cls:=Class(Leader(p,ord),ord); \r
1889    bs:=bs minus {p}; \r
1890         \r
1891    while nops(bs) <>0  do \r
1892      for i in bs do \r
1893         if Class(Leader(i,ord),ord) = cls then poly:=PD_spoly(i,p,ord); res:=res union {poly} fi; \r
1894      od; \r
1895       \r
1896       \r
1897    p:=op(1,bs); \r
1898    cls:=Class(Leader(p,ord),ord); \r
1899    bs:=bs minus {p};   \r
1900       \r
1901    od; \r
1902   \r
1903    res minus {0}; \r
1904 end: \r
1905     \r
1906           \r
1907 \r
1908       \r
1909             \r
1910     \r
1911     \r
1912 PD_spoly:=proc(p,q,ord)     \r
1913 local     leadp,leadq,var,l1,l2,ml,dp,dq,i,Diff_Var; \r
1914 global  depend_relation; \r
1915 options remember,system; \r
1916     leadp:=Leader(p,ord); \r
1917     leadq:=Leader(q,ord); \r
1918      \r
1919     var:=op(0,leadp); \r
1920     Diff_Var:=depend_relation[2]; \r
1921      \r
1922     l1:=dOrder(leadp,var); \r
1923     l2:=dOrder(leadq,var); \r
1924      \r
1925     ml:=maxlist(l1,l2); \r
1926      \r
1927     dp:=p; \r
1928     dq:=q; \r
1929      \r
1930     for i to nops(ml) do \r
1931       dp:=dDiff(dp,Diff_Var[i],ml[i]-l1[i]); \r
1932       dq:=dDiff(dq,Diff_Var[i],ml[i]-l2[i]); \r
1933     od; \r
1934      \r
1935 ######Ö±¾õ¾õµÃ¿ÉÒÔÖ±½ÓÓÃPDPrem(dp,q,Leader(q,ord)) ¸üºÃ¡£¶øÕâÀïÓÃµÄÊÇ±ê×¼µÄ·½·¨¡£      \r
1936    NPremO(dp,dq, Leader(dq,ord)); \r
1937     \r
1938 end:     \r
1939      \r
1940      \r
1941 maxlist:=proc(l1,l2)      \r
1942 local i,list;\r
1943 options remember,system;  \r
1944     list:=[]; \r
1945     for i to nops(l1) do \r
1946       list:=[op(list), max(l1[i],l2[i])]; \r
1947     od: \r
1948 list \r
1949 end:      \r
1950 \r
1951 #torder   Is_deriv Reducedp Leader Head Leader PDE map Num spoly NewPrem \r
1952 \r
1953 # the class of poly f wrt variable ordering ord \r
1954 MainVariables := proc(ps,ord) \r
1955     local i,set; \r
1956        options remember,system; \r
1957     set:={}; \r
1958     for i in ps do \r
1959       set:=set union {Leader(i,ord)}; \r
1960     od; \r
1961 set \r
1962 end:      \r
1963             \r
1964 PNormalp:=proc(p,as,ord) \r
1965     local i,idts,mvs,initp;\r
1966     options remember,system;      \r
1967     initp:=Initial(p,ord); \r
1968     idts:=indets(initp); \r
1969     if idts={} then RETURN(true) fi; \r
1970     mvs:=MainVariables(as,ord); \r
1971     for i in idts do \r
1972        if member(i,mvs) then RETURN(false) fi; \r
1973     od; \r
1974 true \r
1975 end:      \r
1976 \r
1977 ASNormalp:=proc(ps,ord) \r
1978     local i,l,as,p; \r
1979     options remember,system; \r
1980     l:=nops(ps); \r
1981     if l=1 then RETURN(true) fi; \r
1982     as:=[ps[l]]; \r
1983     for i from l-1 to 1 by -1 do \r
1984       p:=ps[i]; \r
1985       if not PNormalp(p,as,ord) then RETURN(false) fi; \r
1986       as:=[p,op(as)]; \r
1987     od;        \r
1988 true \r
1989 end:      \r
1990 \r
1991 Dec:=proc(ps,ord,nzero,T) \r
1992 global remfac,checknum; \r
1993 local i,cset1,cset;\r
1994 options remember,system; \r
1995 \r
1996     cset1:=Cs_c(ps,ord,nzero,T); \r
1997     cset:={}; \r
1998     for i in cset1 do \r
1999       if ASNormalp(i,ord) then \r
2000          cset:=cset union NormDec1(ps,i,ord,nzero,T);      \r
2001          else \r
2002                 cset:=cset union NormDec2(ps,i,ord,nzero,T); \r
2003          fi;                     \r
2004    od; \r
2005 cset \r
2006 end: \r
2007 \r
2008 NormDec1:=proc(ps,as,ord,nzero,T) \r
2009 local     cset,j,nzero1, ps1,ps2,Is,Ninits;\r
2010 options remember,system;\r
2011 \r
2012         cset:= {as};      \r
2013          if Class(as[nops(as)],ord)=1 then Is:=Inits1(as,ord) else Is := Inits(as,ord) fi;             \r
2014          Ninits := Is minus {op(nzero)}; \r
2015      \r
2016          for j to nops(Ninits) do \r
2017                ps1 := ({op(ps)} union {op(as)}) union {Ninits[j]}; \r
2018                if j = 1 then nzero1 := {op(nzero)} \r
2019                else nzero1 := nzero1 union {Ninits[j-1]} \r
2020                fi; \r
2021                cset := cset union Dec(ps1,ord,nzero1,T) \r
2022          od; \r
2023 cset \r
2024 end:           \r
2025           \r
2026 NormDec2:=proc(ps,as,ord,nzero,T) \r
2027 local cset,i,j,l,regas,p,initp,lp,pol,r1,r2,Is,Ninits,l2,ps1,nzero1,mvar,vars,mv,d,f,g;\r
2028 options remember,system; \r
2029     cset:={}; \r
2030      \r
2031     l:=nops(as); \r
2032     regas:=[as[l]]; \r
2033      \r
2034     for i from 2 to l do \r
2035       p:=as[l-i+1]; \r
2036       mvar:=Leader(p,ord); \r
2037       d:=degree(p,mvar); \r
2038       if PNormalp(p,regas,ord) then \r
2039          regas:= [p,op(regas)] \r
2040       else \r
2041          initp:=Initial(p,ord); \r
2042          vars:=indets(initp); \r
2043          lp:=nops(regas); \r
2044          ###1 ¿ªÊ¼111 \r
2045              for j to lp do \r
2046                  pol:=regas[j]; \r
2047                  mv:=Leader(pol,ord); \r
2048                  if member(mv,vars) then \r
2049 #                  ×¢ÒâResultantº¯Êý \r
2050 \r
2051                  r1:=NResultantE(initp,pol, Leader(pol,ord),'f','g'); \r
2052                  r2:=Premas(r1,regas,ord,'std_asc'); \r
2053                  if r2=0 then \r
2054        \r
2055                           if Class(as[nops(regas)],ord)=1 then Is:=Inits1(regas,ord) else Is := Inits(regas,ord) fi;             \r
2056                           Ninits := Is minus {op(nzero)}; \r
2057                \r
2058                        l2:=nops(Ninits); \r
2059                        nzero1:={op(nzero)}; \r
2060 \r
2061                           for j to l2 do \r
2062                               ps1 := ({op(ps)} union {op(as)}) union {Ninits[j]}; \r
2063                               if j <> 1 then \r
2064                                nzero1 := nzero1 union {Ninits[j-1]} \r
2065                               fi; \r
2066                               cset := cset union Dec(ps1,ord,nzero1,T) \r
2067                           od; \r
2068                           print("kkkkkkkkkk"); \r
2069                           if l2 <> 0 then nzero1:=nzero1 union {Ninits[l2]} fi; \r
2070 #                   f;g; \r
2071 #                   print("cset=",cset);\r
2072 #                   print("ps=",ps);\r
2073 #                   print("as=", as);\r
2074 #                   print("regas=",regas);\r
2075                    \r
2076                           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); \r
2077                    RETURN(cset);            \r
2078                  else \r
2079              \r
2080 #                  p ±äÐÎ³ÉÐÂµÄp \r
2081               p:=expand(f*p+g*pol*mvar^d); \r
2082               p:=Premas(p,regas,ord,'std_asc'); \r
2083                  initp:=Initial(p,ord); \r
2084                  vars:=indets(initp); \r
2085               fi; \r
2086               fi; \r
2087              od;   \r
2088          ###1 ½áÊø111   \r
2089          regas:=[p,op(regas)]   \r
2090       fi;        \r
2091       od; \r
2092     cset:=NormDec1(ps,regas,ord,nzero,T); \r
2093 cset      \r
2094 end:        \r
2095 \r
2096 \r
2097 #Èç¹ûÊ×ÏîÏµÊý±ä³É0£¬¶øÓàÊ½²»ÊÇ0ÔòÖÃdegenerateÎª1·ñÔòÎª0¡£\r
2098 WuPrem:=\r
2099 \r
2100    proc(uu,vv,x,ord,degenerate)\r
2101    local r,init,lead,lmono,red,init_r,s,t;\r
2102    options remember,system;    \r
2103    degenerate:=0;\r
2104    if uu=0 then RETURN(0) fi;\r
2105       r:=uu;  \r
2106    if Class(uu,ord) = Class(vv,ord) then \r
2107       r:= NPrem(uu,vv,x);\r
2108    elif Class(uu,ord) > Class(vv,ord) then\r
2109       init:=Initial(uu,ord); \r
2110       if  degree(init,x) >0 and Reducedp(init, vv,ord,'std_asc') = 0  then \r
2111           lead:=Leader(uu,ord);\r
2112 #         collect(uu,lead);\r
2113           lmono:=lead^degree(uu,lead);\r
2114           red:=expand(uu-init*lmono);\r
2115           init_r:=NPremE(init,vv,x,'s','t');    \r
2116           r:=expand(init_r*lmono+s*red); \r
2117           if init_r=0 and r <>0 then degenerate:=1 fi;\r
2118 #          collect(r,lead);\r
2119        fi;\r
2120     fi;\r
2121     r;\r
2122 end:    \r
2123 \r
2124 NPremE := \r
2125 \r
2126    proc(uu,vv,x,s1,t1) \r
2127    local r,v,dr,dv,l,t,lu,lv,s0,t0; \r
2128    options remember,system; \r
2129 \r
2130         s0:=1; \r
2131         t0:=0;           \r
2132            r := expand(uu); \r
2133            dr := degree(r,x); \r
2134            v := expand(vv); \r
2135            dv := degree(v,x); \r
2136            if dv <= dr then l := coeff(v,x,dv); v := expand(v-l*x^dv) \r
2137            else l := 1 \r
2138            fi; \r
2139            while dv <= dr and r <> 0 do \r
2140                gcd(l,coeff(r,x,dr),'lu','lv'); \r
2141          s0:=lu*s0; \r
2142          t0:=lu*t0+lv*x^(dr-dv); \r
2143                t := expand(x^(dr-dv)*v*lv); \r
2144                if dr = 0 then r := 0 else r := subs(x^dr = 0,r) fi; \r
2145                r := expand(lu*r)-t; \r
2146                dr := degree(r,x) \r
2147            od; \r
2148         s1:=expand(s0); \r
2149         t1:=expand(t0);           \r
2150            r \r
2151    end: \r
2152 \r
2153 #This function is for ascending set normalization. \r
2154 #in fact, it is not  a resultant for 2 pols with indeterminate x. \r
2155 #for two polys uu,vv with indeterminate x \r
2156 #it has the property m*uu+n*vv=r with degree(r,x)=0; \r
2157 NResultantE := \r
2158 \r
2159    proc(uu,vv,x,m,n) \r
2160    local r,s,t,r0,s0,t0,r1,s1,t1,tmpr,tmps,tmpt; \r
2161    options remember,system; \r
2162         r0:=uu; \r
2163         s0:=1; \r
2164         t0:=0; \r
2165         r1:=vv; \r
2166         s1:=0; \r
2167         t1:=1;           \r
2168         if degree(r0,x)=0 then m:=1;n:=0;RETURN(r0) fi; \r
2169            if degree(r1,x)=0 then m:=0;n:=1;RETURN(r1) fi; \r
2170 \r
2171            while degree(r1,x) >= 1 do \r
2172              r:=NPremE(r0,r1,x,'s','t');        \r
2173           \r
2174             tmpr:=r1; \r
2175          tmps:=s1; \r
2176          tmpt:=t1; \r
2177          r1:=r; \r
2178          s1:=s*s0-t*s1; \r
2179          t1:=s*t0-t*t1; \r
2180          r0:=tmpr; \r
2181          s0:=tmps; \r
2182          t0:=tmpt; \r
2183            od; \r
2184         m:=s1; \r
2185         n:=t1; \r
2186            r \r
2187    end: \r
2188 \r
2189 \r
2190 \r
2191 ##############################################################################\r
2192 # Name: SylvesterResultant\r
2193 # Calling sequence:  SylvesterResultant(f, g, x, u, v)\r
2194 #   Input: f, g two polynomials in x with positive degrees\r
2195 #          x, a name\r
2196 #   Output: the Sylvester resultant of f and g with respect to x\r
2197 #   Optional: u, v,  such that \r
2198 #                   u*f+v*g = resultant\r
2199 ###############################################################################\r
2200 \r
2201 #SylvesterResultant := proc(f, g, x, u, v)\r
2202 NResultantE := proc(f, g, x, u, v)\r
2203    local df, dg, deg, A, a, b, U, V, i, m, Q, l, e, s, us, vs;\r
2204 \r
2205    #-------------------------------------------\r
2206    # check input\r
2207    #-------------------------------------------\r
2208 \r
2209    (df, dg) := degree(f, x), degree(g, x);\r
2210    if df = 0 or dg = 0 then \r
2211       error "input polynomials should be of positive degrees";\r
2212    end if;\r
2213 \r
2214    #---------------------------------------------\r
2215    # initialize data\r
2216    #---------------------------------------------\r
2217 \r
2218    if df >= dg then\r
2219       (A[1], A[2]) := f, g;\r
2220       (deg[1], deg[2]) := df, dg;\r
2221    else\r
2222       (A[1], A[2]) := g, f;\r
2223       (deg[1], deg[2]) := dg, df;\r
2224    end if;\r
2225 \r
2226    if nargs > 3 then\r
2227          (U[1], U[2]) := 1, 0;\r
2228          if nargs > 4 then\r
2229             (V[1], V[2]) := 0, 1;\r
2230          end if;\r
2231    end if;\r
2232 \r
2233    (a[1], b[1], a[2], b[2]) := (1, 1, lcoeff(A[2],x), lcoeff(A[2], x)^(deg[1]-deg[2]));\r
2234 \r
2235    i := 2;\r
2236    l[i] := deg[i-1]-deg[i]+1;\r
2237 \r
2238    #----------------------------------------------\r
2239    #  compute\r
2240    #----------------------------------------------\r
2241 \r
2242    while true do\r
2243          i := i + 1;\r
2244 \r
2245          #-------------------------------\r
2246          # form extraneous factor\r
2247          #-------------------------------\r
2248 \r
2249          e[i] := (-1)^l[i-1]*a[i-2]*b[i-2];\r
2250 \r
2251          #--------------------------------\r
2252          # compute remainder\r
2253          #--------------------------------\r
2254 \r
2255          if nargs = 3 then\r
2256             A[i] := evala(Simplify(prem(A[i-2], A[i-1], x)/e[i]));\r
2257          else\r
2258             A[i] := evala(Simplify(prem(A[i-2], A[i-1], x, 'm', 'Q')/e[i]));\r
2259             U[i] := evala(Simplify((m*U[i-2] - Q*U[i-1])/e[i]));\r
2260             if nargs = 5 then\r
2261                V[i] := evala(Simplify((m*V[i-2] - Q*V[i-1])/e[i]));\r
2262             end if;\r
2263           end if;\r
2264 \r
2265           #-----------------------------------\r
2266           # Resultant is zero\r
2267           #-----------------------------------\r
2268                     \r
2269           if A[i] = 0 then \r
2270              if nargs > 3 then\r
2271                 u := U[i];\r
2272                 if nargs > 4 then\r
2273                    v := V[i];\r
2274                 end if;\r
2275               end if;\r
2276               return A[i];\r
2277            end if;\r
2278 \r
2279            #------------------------------------\r
2280            # update auxiliary sequences\r
2281            #------------------------------------\r
2282 \r
2283            deg[i] := degree(A[i], x);\r
2284            l[i] := deg[i-1]-deg[i]+1;\r
2285            a[i] := lcoeff(A[i], x);\r
2286            b[i] := evala(Simplify(a[i]^(l[i]-1)/b[i-1]^(l[i]-2)));\r
2287            \r
2288 \r
2289            #-------------------------------------------\r
2290            # Resultant is nonzero\r
2291            #-------------------------------------------\r
2292 \r
2293            if deg[i] = 0 then\r
2294               if nargs > 3 then\r
2295                  s := evala(Simplify(b[i]/a[i]));\r
2296                  us := evala(Simplify(U[i]*s));\r
2297                  if nargs > 4 then\r
2298                     vs  := evala(Simplify(V[i]*s));\r
2299                  end if;\r
2300                  if df >= dg then\r
2301                     u := us;\r
2302                     if nargs > 4 then\r
2303                        v := vs;\r
2304                     end if;\r
2305                  else\r
2306                     u := (-1)^(df*dg)*vs;\r
2307                     if nargs > 4 then\r
2308                        v := (-1)^(df*dg)*us;\r
2309                     end if;\r
2310                  end if;\r
2311               end if;\r
2312               if df >= dg then\r
2313                  return b[i];\r
2314               else\r
2315                  return (-1)^(df*dg)*b[i];\r
2316               end if;\r
2317            end if;\r
2318     end do;\r
2319 end proc:\r
2320    \r
2321 \r
2322 \r
2323 #printf("    Pls type 'with(linalg)' firstly before calling function 'Resultant_E' "); \r
2324 \r
2325 ##The following is to compute the resultant of 2 pols via computing the determinant directly.\r
2326 # input: M is a list of lists , whereby to denote a matrix; \r
2327 #        i,j are two integer numbers; \r
2328 # output:sM is a new list of lists, obtained by removing \r
2329 #        the ith row and jth column of L, \r
2330 SubM:=proc(M,i,j) \r
2331     local row,tmp,sM; \r
2332     if i<1 or i>nops(M) \r
2333          then error "Input row index %1 broke bounds",i;fi; \r
2334     if j<1 or j>nops(M[1]) \r
2335          then error "Input column index %1 broke bounds",j;fi; \r
2336     sM:=[]; \r
2337     tmp:=subsop(i=NULL,M); \r
2338     for row in tmp \r
2339          do row:=subsop(j=NULL,row); \r
2340             sM:=[op(sM),row]; \r
2341          od; \r
2342     sM; \r
2343 end: \r
2344      \r
2345 \r
2346 \r
2347 \r
2348 # with(linalg) firstly \r
2349 # input: A,B are two polynomials, v is a variable; \r
2350 # output:[R,T,U]; R,T,U are three polynomials, where R is the resultant \r
2351 #        of A,B w.r.t v, so that A*T+B*U=R \r
2352 #with(linalg); \r
2353 NResultantE_Matrix:=proc(A,B,v,TA,UB) \r
2354     local m,n,d1,d2,cA,cB,row,i,j,syl,msyl,sM,sign,R,T,U,result; \r
2355 #    #lprint("nargs=",nargs); \r
2356 \r
2357     if (nargs <> 3) and (nargs <> 5) then error "Number of parameters should be 3 or 5" fi; \r
2358      \r
2359     m:=degree(A,v); \r
2360     n:=degree(B,v); \r
2361     if m=0 then error "InPut Polynomial %1 Has No Variable %2",A,v;fi; \r
2362     if n=0 then error "InPut Polynomial %1 Has No Variable %2",B,v;fi; \r
2363   \r
2364 #  to get coefficients of A and B \r
2365     cA:=[]; \r
2366     for d1 from 0 to m \r
2367          do cA:=[coeff(A,v,d1),op(cA)]; \r
2368          od; \r
2369     cB:=[]; \r
2370     for d2 from 0 to n \r
2371          do cB:=[coeff(B,v,d2),op(cB)]; \r
2372          od; \r
2373 \r
2374 #  initialize sylvester matrix \r
2375      \r
2376     syl:=[]; \r
2377     for i from 1 to n         \r
2378          do row:=cA; \r
2379                for j from 1 to i-1 \r
2380               do row:=[0,op(row)];od; \r
2381             for j from 1 to n-i \r
2382                 do row:=[op(row),0];od; \r
2383                syl:=[op(syl),row]; \r
2384             od; \r
2385              \r
2386     for i from 1 to m \r
2387          do row:=cB; \r
2388                for j from 1 to i-1 \r
2389               do row:=[0,op(row)];od; \r
2390             for j from 1 to m-i \r
2391                 do row:=[op(row),0];od; \r
2392                syl:=[op(syl),row]; \r
2393             od; \r
2394      \r
2395     msyl:=linalg[matrix](m+n,m+n,syl); \r
2396      \r
2397 #   to get resultant of A,B w.r.t v \r
2398     R:=linalg[det](msyl); \r
2399     R:=expand(R); \r
2400     if nargs=3 then RETURN(R) fi; \r
2401      \r
2402 #   to get T \r
2403     T:=0; \r
2404     for d1 from 1 to n \r
2405          do \r
2406          sign:=(-1)^(d1+m+n mod 2);           \r
2407          sM:=SubM(syl,d1,m+n); \r
2408          T:=T+(v^(n-d1))*expand(sign*linalg[det](linalg[matrix](m+n-1,m+n-1,sM))); \r
2409          T:=expand(T); \r
2410          od; \r
2411           \r
2412 #  to get U \r
2413         U:=0; \r
2414     for d2 from 1 to m \r
2415          do \r
2416          sign:=(-1)^(d2+m mod 2);      \r
2417          sM:=SubM(syl,d2+n,m+n); \r
2418          U:=U+(v^(m-d2))*expand(sign*linalg[det](linalg[matrix](m+n-1,m+n-1,sM))); \r
2419          U:=expand(U); \r
2420          od; \r
2421         if nargs=5 then TA:=T; UB:=U fi; \r
2422     R; \r
2423 end: \r
2424 \r
2425 ###########################################################################################\r
2426 ###########################################################################################\r
2427 ############## computing partioned-parametric Grobner basis 2005.11.4###################\r
2428 ##ÐèÒª¡¡¡°wsolve¡± ºÍ¡¡¡°project¡±\r
2429 ##ÎÄºó¸¶ÓÐÖ÷Òªº¯Êýµ÷ÓÃµÄÀý×Ó\r
2430 ###########################################################################################\r
2431 ###########################################################################################\r
2432 ###########################################################################################\r
2433 \r
2434 interface(warnlevel=0);\r
2435 with(Groebner): \r
2436 \r
2437 \r
2438 ########pÔÚÔ¼ÊøZero(a11/a12)ÏÂÒ»¶¨²»Îª£°,ordÊÇ²Î±äÔª########## \r
2439 NotZeroUnderCons:=proc(a11,a12,p,ord) \r
2440 option remember; \r
2441     RETURN(IsEmpty({op(a11),p},a12,ord)); \r
2442 end: \r
2443 \r
2444 ########pÔÚÔ¼ÊøZero(a11/a12)ÏÂÒ»¶¨Îª£°,ordÊÇ²Î±äÔª########## \r
2445 IsZeroUnderCons:=proc(a11,a12,p,ord) \r
2446 option remember; \r
2447     RETURN(IsEmpty(a11,{op(a12),p},ord)); \r
2448 end: \r
2449 \r
2450 ########Ô¼Êø¶àÏîÊ½conspolysetÐÎÊ½:[[{a11}£¬{a12}]£¬[{a21}£¬{a22}]]##\r
2451 UnAmbiguous:=proc(conspolyset,var_para::list,term_order) \r
2452 local p,vp,result,a11,a12,a21,a22,newa11,newa22,lc,vars,paras,reduct; \r
2453 \r
2454 vars:=var_para[1]; \r
2455 paras:=var_para[2]; \r
2456 \r
2457 a11:=conspolyset[1][1]; \r
2458 a12:=conspolyset[1][2]; \r
2459 a21:=conspolyset[2][1]; \r
2460 a22:=conspolyset[2][2]; \r
2461 if a22={} then RETURN ({conspolyset}); fi; \r
2462 \r
2463 p:=a22[1]; \r
2464 newa22:=a22 minus {p};\r
2465 \r
2466 if nops(a21) <>0 then p:=numer(normalf(p,a21,term_order(op(vars)))); fi;\r
2467 \r
2468 ##########AAA\r
2469 if (_field_characteristic_ <>0 )        then p :=expand(p) mod _field_characteristic_ fi;\r
2470 \r
2471 if nops(a11) <>0 then p:=numer(normalf(p,a11,term_order(op(paras)))); fi;\r
2472 \r
2473 if (_field_characteristic_ <>0 )        then p :=expand(p) mod _field_characteristic_ fi;\r
2474 \r
2475 \r
2476 #########Èç¹ûpÊÇ0¶àÏîÊ½############ \r
2477 if p=0 then  result:=UnAmbiguous([[a11,a12],[a21,newa22]],var_para,term_order);  RETURN(result); fi; \r
2478 \r
2479 ##########Èç¹ûpÊÇ·Ç0³£Êý##############   \r
2480 if (p <> 0) and ( type(p,'constant') ) then RETURN ( {[conspolyset[1],[{1},{} ]  ]}); fi; \r
2481 \r
2482 #lprint("p="); #lprint(p);  \r
2483 lc:=leadcoeff(p,term_order(op(vars))); \r
2484 #lprint("lc="); #lprint(lc);\r
2485 reduct:=expand(p-lc*leadterm(p,term_order(op(vars)))); \r
2486 vp:=indets(p); \r
2487     \r
2488 #######Èç¹ûpº¬ÓÐ±äÔª############   \r
2489  if (vp intersect {op(vars)}) <> {} then \r
2490      if(type(lc,'constant')=true) then \r
2491         result:=UnAmbiguous([[a11,a12],[a21 union {p}, newa22]],var_para,term_order); \r
2492 \r
2493      elif NotZeroUnderCons(a11,a12,lc,paras) \r
2494          then result:=UnAmbiguous([[a11,a12],[a21 union {p}, newa22]],var_para,term_order); \r
2495              \r
2496      elif IsZeroUnderCons(a11,a12,lc,paras) \r
2497          then result:=UnAmbiguous([[a11,a12],[a21,newa22 union {reduct}]],var_para,term_order); \r
2498            \r
2499      else         \r
2500          lc:=CutNozeroFactors(lc,a12); \r
2501          newa11:=CutPolyByFactor(a11,lc);\r
2502          result:=UnAmbiguous([[{op(newa11),lc},a12],[a21,newa22 union {reduct}]],var_para,term_order) \r
2503                union UnAmbiguous([[a11,a12 union FactorList(lc)],[a21 union {p}, newa22]],var_para,term_order); \r
2504                  \r
2505      fi; \r
2506  else \r
2507 \r
2508 #######Èç¹ûp²»º¬±äÔª£¬Ö»º¬²ÎÊý############   \r
2509 ##modify \r
2510      if NotZeroUnderCons(a11,a12,p,paras) then result:={[conspolyset[1],[{1},{} ]]}; \r
2511      elif IsZeroUnderCons(a11,a12,p,paras) then result:=UnAmbiguous([[a11,a12],[a21,newa22]],var_para,term_order) ;           \r
2512         \r
2513      else \r
2514           lc:=CutNozeroFactors(lc,a12);  \r
2515           newa11:=CutPolyByFactor(a11,lc);    \r
2516           result:=UnAmbiguous([[{op(newa11),lc},a12],[a21,newa22]],var_para,term_order) \r
2517                   union UnAmbiguous([[a11,a12 union FactorList(lc)],[{1}, {}]],var_para,term_order); \r
2518      fi; \r
2519  fi; \r
2520 \r
2521 result: \r
2522 end: \r
2523 \r
2524 \r
2525 \r
2526 \r
2527 \r
2528 \r
2529 ######pol¶Ôlist1ÖÐµÄ¶àÏîÊ½Öð¸ö×ö³ý·¨£¬·µ»ØÉÌ#######\r
2530 DivideList:=proc(pol,list1)\r
2531 local i,q,rt;\r
2532    rt:=pol;\r
2533    for i in list1 do\r
2534       if divide(rt,i,'q') then   rt:=q; fi;\r
2535    od;\r
2536    rt;\r
2537 end:\r
2538 \r
2539 ###<1>list1ÖÐµÄ¶àÏîÊ½ÒÑ²»¿ÉÔ¼\r
2540 ###<2>È¥µôpolÖÐÔÚlist1ÀïÁÐ³öµÄÒò×ÓÇÒºöÂÔÖØÊý\r
2541 CutNozeroFactors:=proc(pol,list1)\r
2542 local i,q,sqf_pol,rt;\r
2543    sqf_pol:=sqrfree(pol);\r
2544    sqf_pol:=sqf_pol[2];\r
2545    for i from 1 to nops(sqf_pol) do\r
2546       sqf_pol[i][1]:=DivideList(sqf_pol[i][1],list1);\r
2547    od;\r
2548    rt:=1;\r
2549    for i from 1 to nops(sqf_pol) do\r
2550       rt:=rt*(sqf_pol[i][1]);\r
2551    od;\r
2552    rt;\r
2553 end:      \r
2554 \r
2555 ###È¥µôlist1ÖÐÒÔpolÎªÒò×ÓµÄ¶àÏîÊ½\r
2556 CutPolyByFactor:=proc(list1,pol)\r
2557 local p,rt,q;\r
2558     rt:={};\r
2559     for p in list1 do\r
2560        if not divide(p,pol,'q') then rt:={op(rt),p};fi;\r
2561     od;\r
2562 rt; \r
2563 end:     \r
2564 \r
2565 \r
2566 \r
2567 ##<1>Í¨¹ýgb¶ÔµÈÓÚÁãµÄÔ¼Êø²¿·Ö×önormalformÀ´»¯¼ò\r
2568 ##<2>ÒòÎªgbµÄÊ×ÏµÊýÒ»¶¨²»ÎªÁã£¬Í¨¹ý¶Ôgb×öprimpartÀ´»¯¼ò\r
2569 SimplifyConsgb:=proc(consgb,var_para,term_order)\r
2570 local rt,i,temp,paras,vars;\r
2571     vars:=var_para[1];\r
2572     paras:=var_para[2];\r
2573     rt:={};\r
2574     for i in consgb do\r
2575         temp:=map(normalf,i[2],i[1][1],term_order(op(paras)));\r
2576         temp:=map(primpart,map(numer,temp),vars);        \r
2577         rt:={op(rt),[i[1],temp]};\r
2578     od;\r
2579     rt;\r
2580 end:\r
2581 \r
2582 \r
2583 FactorList:=proc(pol) \r
2584 local i,list1,list2; \r
2585   list1 := factors(pol)[2]; \r
2586   list2 := {}; \r
2587   for i in list1 do  list2 := {op(list2),i[1]}; od; \r
2588   list2 \r
2589 end: \r
2590  \r
2591 \r
2592 #######Èç¹ûpÔÚconstrainµÄÌõ¼þÏÂÎª0ÔÚ·µ»Ø1,np=0 \r
2593 #######Èç¹ûpÔÚconstrainµÄÌõ¼þÏÂÎª·Ç0µÄ³£Êý·µ»Ø10,np=10  \r
2594 #######Èç¹ûpÔÚconstrainµÄÌõ¼þÏÂÊ×ÏîÏµÊýÒ»¶¨²»ÊÇ0Ôò·µ»Ø-1 \r
2595 #######ÆäËûÇéÐÎ·µ»Ø0¡£np=pol \r
2596 ########paracons=[ [part=0 ],[part <>0]] ÊÇ²ÎÊýµÄÔ¼Êø; \r
2597 \r
2598 nonzerolc:=proc(pol,paracons,var_para,np,term_order) \r
2599 local lc,vars,paras,cons1,cons2,np0,t,pol2; \r
2600 global _field_characteristic_ ;\r
2601     vars:=var_para[1]; \r
2602     paras:=var_para[2]; \r
2603     cons1:=paracons[1]; \r
2604     cons2:=paracons[2];\r
2605     lc:=leadcoeff(pol,term_order(op(vars))); \r
2606     if pol=0 then np:=0; RETURN(1); fi; \r
2607     if type(pol,'constant') then np:=10; RETURN(10); fi;       \r
2608     if type(lc,'constant') then np:=pol; RETURN(-1); fi; \r
2609      \r
2610     \r
2611     ##########lcÔÚÔ¼ÊøÏÂÒ»¶¨Îª£°#######   \r
2612     if IsZeroUnderCons(cons1,cons2,lc,paras) then \r
2613          pol2:=reductum(pol,term_order(op(vars))); \r
2614          if (_field_characteristic_ <>0  ) then pol2:=pol2 mod _field_characteristic_ fi;\r
2615          t:=nonzerolc(pol2,paracons,var_para,'np0',term_order); \r
2616          np:=np0; \r
2617          RETURN(t); \r
2618           \r
2619    ##########lcÔÚÔ¼ÊøÏÂÒ»¶¨²»Îª£°####### \r
2620     elif NotZeroUnderCons(cons1,cons2,lc,paras) then  \r
2621          if lc=pol then np:=10; RETURN(10);\r
2622          else  np:=pol;  RETURN(-1);\r
2623          fi; \r
2624     \r
2625     else \r
2626         np:=pol; \r
2627         RETURN(0); \r
2628     fi; \r
2629 end: \r
2630 \r
2631 reductum:=proc(p,ord) \r
2632 local res; \r
2633     expand(p-leadcoeff(p,ord)*leadterm(p,ord)); \r
2634 end:       \r
2635 \r
2636 gb:=proc(ps,var_para,term_order) \r
2637 local res; \r
2638   res:=gb0([[{},{}],[{},{op(ps)}]],var_para,term_order); \r
2639   res:=SimplifyConsgb(res,var_para,term_order);\r
2640 end:   \r
2641 \r
2642 \r
2643 \r
2644 gb0:=proc(conspolyset,var_para,term_order) \r
2645 local i,branches, res,temp; \r
2646    res:={}; \r
2647    branches:=UnAmbiguous(conspolyset,var_para,term_order); \r
2648    for i in branches do \r
2649         temp:=gb1(i,var_para,term_order);\r
2650         res:=res union temp; \r
2651    od; \r
2652    res; \r
2653 end:   \r
2654 \r
2655 #################Step 1######################## \r
2656 gb1:=proc(conspolyset,var_para,term_order) \r
2657 local i,j,T,paracons,pset,pset0,pset1,pol,pset3,l,nonzerosp,sp,bs2,vars,newpol,cr,nf,np; \r
2658 global _field_characteristic_;\r
2659         T:={}; \r
2660         vars:=var_para[1]; \r
2661         paracons:=conspolyset[1];\r
2662         pset:=conspolyset[2][1]; \r
2663           \r
2664         pset0:=pset; \r
2665         pset1:={}; \r
2666           \r
2667         if nops(pset)=0 then  RETURN({[paracons,[0]]}) fi; \r
2668           \r
2669         if  nops(pset) =1 then     RETURN({[paracons,[op(pset)]]}) fi; \r
2670         \r
2671         #################Step 2-----  Normalization------#####################             \r
2672         for i in pset do \r
2673                pset0:=pset0 minus {i}; \r
2674                nf:=normalf(i,pset0 union pset1,term_order(op(vars))); \r
2675                pol:=expand(numer(nf )); \r
2676                cr:=nonzerolc(pol,paracons,var_para,'newpol',term_order); \r
2677                \r
2678                newpol:=numer(normalf(newpol,pset0 union pset1,term_order(op(vars))));\r
2679 \r
2680                if (_field_characteristic_ <>0 ) then newpol :=expand(newpol) mod _field_characteristic_ fi;\r
2681                #lprint("newpol=");#lprint(newpol);\r
2682                if cr=1 then ; \r
2683                elif cr=10 then RETURN( { [paracons,[1]] }); \r
2684                elif cr=-1 then if   newpol <>0 then pset1:=pset1 union {newpol} fi; \r
2685                else \r
2686                   pset3:=[paracons,[ pset0 union pset1 ,{newpol} ]]; \r
2687                   RETURN( gb0(pset3,var_para,term_order) ); \r
2688                 fi;                       \r
2689         od; \r
2690         \r
2691         ##############Step 3-----  S-polynomial------#################################### \r
2692         l:=nops(pset1); \r
2693         if  l=1 then RETURN( {[paracons,[op(pset1)]] }) fi; \r
2694         \r
2695         nonzerosp:={}; \r
2696         for i to l do \r
2697             for j from i+1 to l do \r
2698                sp:=expand(numer(normalf(spoly(pset1[i],pset1[j],term_order(op(vars))),pset1,term_order(op(vars))))); \r
2699                if( _field_characteristic_ <>0) then sp:= expand(sp) mod _field_characteristic_ fi;\r
2700                       #lprint("sp="); #lprint(sp);\r
2701                if (nonzerolc(sp,paracons,var_para,'np',term_order)<>1) then  nonzerosp:=nonzerosp union {np} fi; \r
2702             od; \r
2703         od; \r
2704         \r
2705         if nonzerosp={} then RETURN( {[paracons,inter_reduce([op(pset1)],term_order(op(vars)) )] } ) ; \r
2706         else \r
2707                  bs2:=[paracons,[{op(pset1)},nonzerosp ]]; \r
2708                  T:=T union gb0(bs2,var_para,term_order); \r
2709         fi; \r
2710 \r
2711         T: \r
2712 end:\r
2713 \r
2714 \r
2715 #############################\r
2716 ###   solution number    ####\r
2717 #############################\r
2718 \r
2719 basis:=proc(a,b,ord)\r
2720 local i,j,k1,k;k:={};\r
2721 for i in a do\r
2722    k1:=normalf(i,b,ord);\r
2723    if(k1<>0) then \r
2724      k:=k union {k1};\r
2725    end if;  \r
2726 end do;\r
2727 return k;\r
2728 end:   \r
2729 \r
2730 mulset:=proc(a,b)\r
2731 local i,j,k;\r
2732 k:={};\r
2733 if b={} then return a; fi;\r
2734 for i in a do\r
2735   for j in b do\r
2736      k:=k union {i*j};\r
2737   od;\r
2738 od;\r
2739 return k;\r
2740 end:    \r
2741 \r
2742 create:=proc(a,b,num)\r
2743 local i,j,k,k1,b1;\r
2744 \r
2745 b1:=b union {1}; \r
2746 k1:=b1 union mulset(a,b1);\r
2747 if num=0 then return k1;fi;\r
2748   k:=create(a,k1,num-1);\r
2749 \r
2750 end:    \r
2751 \r
2752 havefinitesolution:=proc(l,var)\r
2753 local i,lvs,lv1;\r
2754 lv1:={};\r
2755 for i in l do\r
2756    lvs:=indets(i) intersect {op(var)};\r
2757    if nops(lvs)=1 then lv1:=lv1 union lvs; fi;\r
2758 od;\r
2759 \r
2760 if ({op(var)} minus lv1 )={}  then return true;\r
2761 else return false;\r
2762 fi;\r
2763 end: \r
2764  \r
2765   \r
2766 comp:=proc(a,b)\r
2767 local i,j,k,l;k:={};\r
2768 for i in b do\r
2769   for j in a do\r
2770      if j=i then   \r
2771         k:=k union {j};\r
2772      fi;\r
2773    od;     \r
2774 od;\r
2775 return k;\r
2776 end:   \r
2777 \r
2778 getmostdegree:=proc(l)\r
2779 local i,j,k1,k2;\r
2780 k2:={};\r
2781 k1:={};\r
2782 for i in l do \r
2783   k1:={op(i)};\r
2784   for j in k1 do\r
2785     k2:=k2 union {op(j)};\r
2786   od;\r
2787 od;    \r
2788 return k2;\r
2789 end:\r
2790 \r
2791 ltt:=proc(l,vorder)\r
2792  local i,k;\r
2793  k:={};\r
2794  for i to nops(l) do\r
2795  k:={leadterm(l[i],vorder)} union k;\r
2796  end do;\r
2797  return k;\r
2798  end:\r
2799 \r
2800 GetDim:=proc(lts,vars)\r
2801 local i,branch,l,rt;\r
2802 rt:=nops(lts);\r
2803 branch:=Nrs(lts,[op(vars)],{});\r
2804 branch:=[op(branch)];\r
2805 for i to nops(branch) do\r
2806    l:=nops(branch[i]);\r
2807    if l<rt then rt:=l;fi;\r
2808  od;\r
2809  return nops(vars)-rt;\r
2810  end:\r
2811 \r
2812 \r
2813 ##for reduced groebner basis\r
2814 sumofdegree:=proc(lts)\r
2815 local i,rt;\r
2816 rt:=0;\r
2817 for i in lts do\r
2818    if nops(indets(i))=1 then\r
2819        rt:=rt+degree(i);\r
2820    fi;\r
2821 od;\r
2822 return rt;\r
2823 end: \r
2824 \r
2825 solution1:=proc(ps,ord,term_order)\r
2826 local i,j,k,cgb,lt,mosttotaldeg,dim;\r
2827 cgb:=gb(ps,ord,term_order);\r
2828 cgb:=[op(cgb)];\r
2829 for i to nops(cgb) do\r
2830   lt:=ltt(cgb[i][2],term_order(op(ord[1])));\r
2831   if  havefinitesolution(lt,ord[1]) then\r
2832     mosttotaldeg:=sumofdegree(lt); \r
2833     k:=basis(create({op(ord[1])},{},mosttotaldeg),lt,term_order(op(ord[1])));\r
2834     cgb[i]:=[op(cgb[i]),k,cat("finite solution with number: ", nops(k))];\r
2835 \r
2836   elif cgb[i][2]=[1] then cgb[i]:=[op(cgb[i]),"no solution"];\r
2837   elif cgb[i][2]=[0] then cgb[i]:=[op(cgb[i]),cat("infinite solution with dimension: ",nops(ord[1]))];\r
2838   else \r
2839        dim:=GetDim(lt,ord[1]);\r
2840        cgb[i]:=[op(cgb[i]),[lt],cat("infinite solution with dimension: ", dim)];\r
2841   fi;\r
2842 od;\r
2843 return {op(cgb)};\r
2844 end:\r
2845 \r
2846 \r
2847 ##µÚ4¸ö²ÎÊý¿ÉÊ¡,Èç¹ûÓÐ,±ØÐë¼Óµ¥ÒýºÅ!\r
2848 solution:=proc(ps,ord,term_order,pgb)\r
2849 local i,j,k,cgb,lt,mosttotaldeg,dim;\r
2850 cgb:=gb(ps,ord,term_order);\r
2851 cgb:=[op(cgb)];\r
2852 if nargs=4 then pgb:={op(cgb)};fi;\r
2853 for i to nops(cgb) do\r
2854   lt:=ltt(cgb[i][2],term_order(op(ord[1])));\r
2855   if  havefinitesolution(lt,ord[1]) then\r
2856     mosttotaldeg:=sumofdegree(lt); \r
2857     k:=basis(create({op(ord[1])},{},mosttotaldeg),lt,term_order(op(ord[1])));\r
2858     cgb[i][2]:=cat("finite solution with number: ", nops(k)); \r
2859   elif cgb[i][2]=[1] then cgb[i][2]:="no solution";\r
2860   elif cgb[i][2]=[0] then cgb[i][2]:=cat("infinite solution with dimension: ",nops(ord[1]));\r
2861   else \r
2862        dim:=GetDim(lt,ord[1]);\r
2863        cgb[i][2]:= cat("infinite solution with dimension: ", dim);\r
2864   fi;\r
2865 od;\r
2866 {op(cgb)};\r
2867 end:\r
2868 \r
2869 FindRightGB:=proc(paragb)\r
2870 local gb,pgb,rt,i,j;\r
2871 j:=0;\r
2872 pgb:=[op(paragb)];\r
2873 for i to nops(pgb) do\r
2874    gb:=pgb[i];\r
2875    if nops(gb[1][1])=0 then gb[1][1]:={0} fi;\r
2876     if nops(gb[1][2])=0 then gb[1][2]:={1} fi;\r
2877    if {op(gb[1][1])}={0} and (gb[1][2] intersect {0})<>{0} then \r
2878            j:=j+1;\r
2879            rt:=gb[2];\r
2880    fi;\r
2881 od;\r
2882 if j<>1 then print("ambiguous occur!");fi;\r
2883 rt;\r
2884 end:\r
2885 \r
2886 EvalAll:=proc(op,vars,vals)\r
2887 local i,rt;\r
2888 rt:=op;\r
2889 for i to nops(vars) do\r
2890     rt:=eval(rt,vars[i]=vals[i]);\r
2891 od;\r
2892 rt;\r
2893 end:\r
2894 \r
2895 \r
2896 #°ÑpsÖÐµÄ¶àÏîÊ½Éè³ÉÊ×ÏîÏµÊýÎª£±\r
2897 SetLC1:=proc(ps)\r
2898 local i,p,rt;\r
2899 rt:={};\r
2900 for i to nops(ps) do \r
2901     p:=expand(ps[i]);\r
2902     if p=0 then rt:={op(rt),p};\r
2903     else\r
2904          p:=p/lcoeff(p);   \r
2905          rt:={op(rt),p};\r
2906     fi;\r
2907 od;\r
2908 rt;\r
2909 end:\r
2910 \r
2911 \r
2912 ###¼ìÑéÔ¼Êøgroebnerbasis³ÌÐòµÄÕýÈ·ÐÔ  \r
2913 check:=proc(paravalue,paragb,ps,ord,term_order)\r
2914 local i,gb,gb1,spec_paragb,spec_ps;\r
2915 if nops(paravalue)<>nops(ord[2]) then return "Error in parameter assignment: parameter number is NOT matched!" fi;   \r
2916 \r
2917 spec_ps:=EvalAll(ps,ord[2],paravalue);\r
2918 spec_ps:={op(spec_ps)} minus {0};\r
2919  spec_ps:=[op(spec_ps)];\r
2920 gb1:=gbasis(spec_ps,term_order(op(ord[1])));\r
2921 \r
2922 spec_paragb:=EvalAll(paragb,ord[2],paravalue);\r
2923 gb:=FindRightGB(spec_paragb);\r
2924 \r
2925 if SetLC1(gb)=SetLC1(gb1) then print( "two results are the same polynomial set!");\r
2926 else print("two results are NOT the same polynomial set!");\r
2927 fi;\r
2928 \r
2929 print("greobner basis gotten by parameter assignment at begin:");\r
2930 print(gb1);\r
2931 \r
2932 print("greobner basis gotten by CGB method:");\r
2933 gb;\r
2934 \r
2935 end:\r
2936 \r
2937 \r
2938 ########  Àý×Ó ######################\r
2939 #²ÎÊý£º\r
2940 #ps:=[a*x^2-b*x,b*y^2-c*x,c*z^2-a*x+b];\r
2941 #ord:=[[x,y,z],[a,b,c]];\r
2942 #termord:=tdeg;\r
2943 #'pgb';¡¡ÊÇÒ»¸ö±äÁ¿Ãû£¬¸Ã²ÎÊý¿ÉÈ±Ê¡\r
2944 \r
2945 #µ÷ÓÃ£±£º\r
2946 #gb(ps,ord,termord)\r
2947 #·µ»Ø£º\r
2948 #psµÄÔ¼ÊøGroebnerBasis\r
2949 \r
2950 #µ÷ÓÃ£²£º\r
2951 #solution(ps,ord,termord,'pgb');¡¡\r
2952 #·µ»Ø£º\r
2953 #psµÄ²»Í¬Ô¼ÊøÏÂ½âµÄ¸öÊý¡£\r
2954 #pgbÎªpsµÄÔ¼ÊøGroebnerBasis\r
2955 \r
2956 #µ÷ÓÃ£³£º\r
2957 #solution1(ps,ord,termord);¡¡\r
2958 #·µ»Ø£º\r
2959 #psµÄÔ¼ÊøGroebnerBasis¼°ÏàÓ¦½âµÄ¸öÊý£¨ÈôÓÐÏÞ£¬»¹ÁÐ³ö²»ÊôÓÚÊ×ÀíÏëµÄµ¥ÏîÊ½£©¡£\r
2960 \r
2961 #²ÎÊý£º\r
2962 #paraval:=[1,2,3];¡¡ÊÇ²ÎÊý[a,b,c]µÄÌØ¶¨Öµ\r
2963 #pgb; psµÄÔ¼ÊøGroebnerBasis\r
2964 \r
2965 #µ÷ÓÃ4£º\r
2966 # check(paraval,pgb,ps,ord,termord); \r
2967 #¸Ãº¯Êý±È½ÏpgbÔÚÌØ¶¨²ÎÊýÖµparavalÏÂµÄGroebnerBasisÓëÏÈ¶Ô²ÎÊýÈ¡ÖµºóÔÙÇóµÄGroebnerBasisÊÇ·ñÏàÍ¬¡£\r
2968 #Í¬Ê±Ò²ÁÐ³öÕâÁ½¸öGroebnerBasis¡£´Ëº¯Êý¿É×öÎªÈ«ÎÄµÄ²âÊÔº¯Êý¡£\r
2969 #·µ»Ø£º\r
2970 #pgbÔÚÌØ¶¨²ÎÊýÖµparavalÏÂµÄGroebnerBasis\r
2971 \r
2972 \r
2973 ###########################################################################################\r
2974 ###########################################################################################\r
2975 ############## computing partioned-parametric Grobner basis 2005.11.4 #####################\r
2976 ############## END END END END END END END END END END END END END END#####################\r
2977 ###########################################################################################\r
2978 ###########################################################################################\r
2979 ###########################################################################################\r
2980 \r
2981 \r
2982 ########################for computing the projection of quasi-variety 2005.11.4###############\r
2983 ###  read "wsolve.txt" firstly  ### \r
2984 if not type(wsolve,procedure) then \r
2985 #lprint(" Warning: please read 'wsolve.txt' firstly ! "):\r
2986 fi:\r
2987 \r
2988 ##-----main function------------------------------------------------\r
2989 ##Project(ps,ds,ord,pord) //old main function\r
2990 ##Sproject(result_of_project,pord) //no component is redundant \r
2991 ##Example:  [[ASC,[a,b]],...] == (ASC=0 and a*b<>0 ) or ...\r
2992 \r
2993 ##Proj(ps,ds,ord,pord)  //new main function \r
2994 ##Sproj(result_of_project,pord)  //no component is redundant \r
2995 ##Example: [[ASC,[a,b]],...] == (ASC=0 and ( a<>0 or b<>0) ) or ...\r
2996 ##-------------------------------------------------------------------\r
2997 \r
2998 \r
2999 #-----------Part1------<algebraic case>---------------\r
3000 #         project(ps,ds,ord,pord)\r
3001 #         ord:= all variables (descending order list)\r
3002 #         pord:=variables to eliminate (as above)\r
3003 #-----------------------------------------------------\r
3004 #QV (QuasiVariety) ; ASC (ascendant characteristic set)\r
3005 \r
3006 # get the product of the initials of cs\r
3007 \r
3008 Mainvar := proc(f,ord)\r
3009                   local i;\r
3010                   options remember,system;\r
3011                       for i to nops(ord) do\r
3012                           if has(f,ord[i]) then RETURN(ord[i]) fi\r
3013                       od;\r
3014                       #lprint(`Warning: lvar is called with constant`);\r
3015                       0\r
3016                   end:\r
3017 \r
3018 GetIP:=proc(cs,ord)\r
3019         local i,rt;\r
3020         rt:=1;\r
3021         for i from 1 to nops(cs) do\r
3022                 rt:=rt*Initial(cs[i],ord);\r
3023         od;\r
3024         rt;\r
3025 end:\r
3026 \r
3027 # get the product of the polys in ds.\r
3028 GetProduct:=proc(ds)\r
3029         local jds,d;\r
3030         jds:=1;\r
3031         for d in ds do\r
3032                 jds:=jds*d;\r
3033                 od;\r
3034         jds;\r
3035 end:\r
3036 \r
3037 #get the union of two lists, denoted as list also.\r
3038 ListUnion:=proc(list1,list2)\r
3039         local rt;\r
3040         rt:={op(list1)} union {op(list2)};\r
3041         rt:=[op(rt)];\r
3042 end:\r
3043 \r
3044 # inverse a list\r
3045 Inverse:=proc(lis)\r
3046     local i,result;\r
3047     result:=[];\r
3048     for i from nops(lis)  to 1 by -1\r
3049         do result:=[op(result),lis[i]];\r
3050         od;\r
3051     result;\r
3052 end:     \r
3053 \r
3054 # inverse every list in a list.\r
3055 InverseAll:=proc(lislis)\r
3056         local rt,l;\r
3057         rt:=[];\r
3058         for l in lislis do\r
3059                 rt:=[op(rt),Inverse(l)];\r
3060         od;\r
3061 end:\r
3062 \r
3063 # eliminate such ASC in cslist that Premas(jds,ASC,ord)=0\r
3064 Newcslist:=proc(cslist,jds,ord)\r
3065         local rt,cs;\r
3066         rt:=[];\r
3067         for cs in cslist do\r
3068                 if Premas(jds,cs,ord)<>0\r
3069                         then rt:=[op(rt),cs]                               \r
3070                 fi;\r
3071                 od;\r
3072         rt;\r
3073 end:\r
3074 \r
3075 # factors of a polynomial w.r.t ord; \r
3076 # facts(x^2*y+x*y-x,[x,y]); [x,x*y+y-1]\r
3077 # facts(x^2*y+x*y-x,[y]); [x*y+y-1]\r
3078 # facts(x^2*y+x*y-x,[u,v,w]); []\r
3079 Ordfacts:=proc(pol,ord)\r
3080    local i,temp,facl,result;\r
3081    if type(pol,polynom) then      \r
3082            result:=[];\r
3083            facl:=factors(pol);\r
3084            facl:=facl[2];\r
3085            if facl=[] then RETURN([]) fi;\r
3086            for i from 1 to nops(facl)\r
3087                do temp:=facl[i];\r
3088                   if dClass(temp[1],ord)<>0 then\r
3089                   result:=[op(result),temp[1]];\r
3090                   fi;\r
3091                od;\r
3092            RETURN(result);\r
3093    else  \r
3094          #lprint(`Warning: facts is called with no polynomial`);\r
3095          RETURN(pol);\r
3096    fi;  \r
3097 end:  \r
3098 \r
3099 # put p to proper position in notzerolist s.t. MV(ASC[i-1])<=MV(p)<MV(ASC[i])\r
3100 PutToNotZero:=proc(p,notzerolist,mvl,ord)\r
3101         local i,cp,rt;\r
3102         rt:=notzerolist;\r
3103         cp:=Class(p,ord);       \r
3104         if cp=0 then RETURN(rt);fi;\r
3105         for i from 1 to nops(mvl) do \r
3106                 if cp<mvl[i] then rt:=subsop(i=rt[i]*p,rt);break;\r
3107                 fi; \r
3108         od;\r
3109         i;\r
3110         if i=nops(mvl)+1 then \r
3111                 rt:=subsop(-1=rt[-1]*p,rt) ;\r
3112         fi;\r
3113         rt;\r
3114 end:    \r
3115 \r
3116 # all factors(noted fpolys) of polys in ds and initials, \r
3117 # notzerolist[i]= the product of p, \r
3118 # s.t. <1> p in fpolys and <2> MV(ASC[i-1])<=MV(p)<MV(ASC[i])       \r
3119 CreatNotZeroList:=proc(cs,ds,mvl,ord)\r
3120         local i,fs,p,notzerolist;\r
3121         \r
3122         notzerolist:=[];\r
3123         for i from 0 to nops(cs) do\r
3124              notzerolist:=[op(notzerolist),1];\r
3125         od;\r
3126         \r
3127         fs:=[];\r
3128         for p in ds do\r
3129                 fs:=ListUnion(fs,Ordfacts(p,ord));\r
3130                 od;\r
3131         for p in fs do\r
3132                 notzerolist:=PutToNotZero(p,notzerolist,mvl,ord);\r
3133                 od;\r
3134         notzerolist;\r
3135 end:\r
3136 \r
3137 #a front part of ASC whose polynomials contain no variables in pord     \r
3138 SubASC:=proc(asc,pord)\r
3139    local i, result;\r
3140    result:=[];\r
3141    for i from 1 to nops(asc)\r
3142        do if Class(asc[i],pord)=0 \r
3143              then result:=[op(result),asc[i]];\r
3144           else break;          \r
3145           fi;\r
3146        od;\r
3147    RETURN(result);\r
3148 end:     \r
3149 \r
3150 #judge the result of ProjectASC\r
3151 #if it is [[...],1] return 1;\r
3152 #if it is [[...],0] return 0;\r
3153 #else return 3;\r
3154 IsTrival:=proc(qv,ord)\r
3155         local rt;\r
3156         if nops(qv)<>2 then rt:=3;\r
3157         elif not type(qv[2],polynom) then rt:=3;\r
3158         elif Class(qv[2],ord)<>0 then rt:=3;\r
3159         elif qv[2]=0 then rt:=0;\r
3160         else rt:=1;\r
3161         fi;  \r
3162         RETURN(rt);\r
3163 end:\r
3164 \r
3165 # get variables of pol which are of higher order then the highest\r
3166 # polynomial in asc. \r
3167 HigherOrdVars:=proc(pol,mvl,len,ord)\r
3168         local vl,rt,var,index;\r
3169         vl:=indets(pol) intersect {op(pord)};\r
3170         rt:=[];\r
3171         if len<=1 then RETURN(rt);fi;\r
3172         for var in vl do\r
3173                 index:=Class(var,ord);\r
3174                 if index>mvl[len-1] then\r
3175                         rt:=[op(rt),var];\r
3176                 fi;\r
3177         od;\r
3178         rt;\r
3179 end:\r
3180 \r
3181 # decompose the pol in [[ASC],pol] to ['UN',[[ASC],pol1],...]\r
3182 # so as to eliminate variables in varl.  \r
3183 DecompList:=proc(qv,subasc,notzerolist,varList,ord)\r
3184         local i,flag,lenasc,rt,pol,coeffsList,coeff;\r
3185         if varList=[] then RETURN(qv);fi;\r
3186         pol:=qv[2];\r
3187         \r
3188         if pol=0 then rt:=[subasc,pol]; \r
3189         else \r
3190                 flag:=0; \r
3191                 lenasc:=nops(qv[1]);\r
3192                 for i from 1 to lenasc do\r
3193                         if notzerolist[i]!=1 then flag:=1; break;fi;\r
3194                 od;\r
3195                 coeffsList:=[coeffs(expand(pol),varList)];\r
3196                 if nops(coeffsList)>1 then rt:=["UN"];\r
3197                 else rt:=[]; \r
3198                 fi;\r
3199                 for coeff in coeffsList do\r
3200                         if Class(coeff,ord)=0 and flag=0 then \r
3201                                 RETURN([subasc,1]);                             \r
3202                         else\r
3203                                 rt:=[op(rt),[qv[1],coeff]];\r
3204                         fi;\r
3205                 od;             \r
3206         fi;\r
3207         if nops(rt)=1 then rt:=rt[1];fi;\r
3208         rt;\r
3209 end:\r
3210 \r
3211 \r
3212 #when ASC has no variables in varl . \r
3213 # decompose the pol in [[ASC],pol] to ['UN',[[ASC],pol1],...]\r
3214 # so as to eliminate variables in varl \r
3215 DecompList1:=proc(qv,subasc,notzerolist,varList,ord)\r
3216         local i,flag,lenasc,rt,pol,coeffsList,coeff;\r
3217         if varList=[] then RETURN(qv);fi;\r
3218         pol:=qv[2];\r
3219         \r
3220         if pol=0 then rt:=[subasc,pol]; \r
3221         else \r
3222                 flag:=1; \r
3223                 lenasc:=nops(qv[1]);\r
3224                 for i from 1 to lenasc do\r
3225                         flag:=flag*notzerolist[i] ; \r
3226                 od;\r
3227                 \r
3228                 coeffsList:=[coeffs(expand(pol),varList)];\r
3229                 if nops(coeffsList)>1 then rt:=["UN"];\r
3230                 else rt:=[]; \r
3231                 fi;\r
3232                 for coeff in coeffsList do\r
3233                         if Class(coeff,ord)=0 and flag=1 then \r
3234                                 RETURN([subasc,1]);                             \r
3235                         else\r
3236                                 rt:=[op(rt),[qv[1],coeff*flag]];\r
3237                         fi;\r
3238                 od;             \r
3239         fi;\r
3240         if nops(rt)=1 then rt:=rt[1];fi;\r
3241         rt;\r
3242 end:\r
3243 \r
3244 \r
3245 #  eliminate the highest polynomial in ASC, within [[ASC],pol] \r
3246 ElimPoly:=proc(qv,notzerolist,mvl,subasc,ord,pord)\r
3247         local rt,newasc,asc,pol,hPol,hVar,hDeg,rem,newPol,tempList,varList,len;\r
3248         asc:=qv[1];\r
3249         len:=nops(asc);\r
3250         newasc:=subsop(len=NULL,asc);\r
3251         pol:=qv[2];\r
3252         hPol:=asc[len];\r
3253         hVar:=Mainvar(hPol,ord);\r
3254         hDeg:=degree(hPol,hVar);\r
3255         if degree(pol,hVar)=0 then \r
3256                 rt:=[newasc,pol*notzerolist[len]];\r
3257                 RETURN(rt);\r
3258         else\r
3259                 rem:=sprem(pol^hDeg,hPol,hVar);\r
3260                 if rem=0 then \r
3261                         rt:=[subasc,0];\r
3262                         RETURN(rt);                     \r
3263                 fi; \r
3264                 newPol:=rem*notzerolist[len];\r
3265                 tempList:=[newasc,newPol];\r
3266                 varList:=HigherOrdVars(newPol,mvl,len,ord,pord);        \r
3267                 rt:=DecompList(tempList,subasc,notzerolist,varList,ord);\r
3268                 RETURN(rt);\r
3269         fi;\r
3270 end:\r
3271 \r
3272 # project wth one QV formed by one ASC(prepair share done by ProjASC)\r
3273 ProjectASC:=proc(qv,notzerolist,mvl,subasc,lensub,ord,pord)\r
3274         local pol,asc,nproj,tempList,sign,rt,nozero,i;\r
3275         pol:=qv[2];\r
3276         asc:=qv[1];\r
3277         \r
3278 #  get the nonzero polynomial to prject with here.    \r
3279         nozero:=pol;\r
3280         for i from 1 to nops(asc) do\r
3281                 nozero:=nozero*notzerolist[i];\r
3282                 od;\r
3283 \r
3284 #  recursion ending conditions          \r
3285         if nops(asc)=0 or nops(asc)=lensub then\r
3286                 if Class(nozero,pord)>0 then \r
3287                         rt:=DecompList1(qv,subasc,notzerolist,pord,ord);\r
3288                 elif Class(nozero,ord)>0 then \r
3289                         rt:=[asc,nozero];\r
3290                 else \r
3291                         rt:=[asc,Get1or0(nozero)];\r
3292                 fi;\r
3293                 RETURN(rt);\r
3294         elif Class(nozero,pord)=0 then\r
3295                 if Class(nozero,ord)>0 then rt:=[subasc,nozero];\r
3296                 else \r
3297                         rt:=[subasc,Get1or0(nozero)];\r
3298                 fi;\r
3299                 RETURN(rt);\r
3300         fi;\r
3301         \r
3302 #  eliminate the highest polynomial in asc of QV.\r
3303         tempList:=ElimPoly(qv,notzerolist,mvl,subasc,ord,pord);\r
3304         \r
3305         # whether have many branchs in the result of ElimPoly.\r
3306         if not type(tempList[1],string) then \r
3307                 rt:=ProjectASC(tempList,notzerolist,mvl,subasc,lensub,ord,pord);\r
3308         else \r
3309                 rt:=[tempList[1]];\r
3310                 for i from 2 to nops(tempList) do\r
3311                         nproj:=ProjectASC(tempList[i],notzerolist,mvl,subasc,lensub,ord,pord);\r
3312                         \r
3313                         # deal with trival cases.\r
3314                         sign:=IsTrival(nproj,ord);\r
3315                         if sign=0 then next;\r
3316                         elif sign=1 then rt:= [subasc,1];break;\r
3317                         else rt:=[op(rt),nproj];\r
3318                         fi;\r
3319                 od:\r
3320                 \r
3321                 # rt:=["UN"]\r
3322                 if nops(rt)=1 then rt:=[subasc,0]; fi;\r
3323                 \r
3324                 # rt:=["UN",[...]]\r
3325                 if nops(rt)=2 and type(rt[2],list) then rt:=rt[2];fi;\r
3326         fi;\r
3327         rt;\r
3328 end: \r
3329 \r
3330 # project with one QV formed by one ASC\r
3331 ProjASC:=proc(cs,ds,ord,pord)\r
3332         local mvl,notzerolist,p,qv,subasc,lensub,rt;\r
3333         \r
3334         # mainvars order list\r
3335         mvl:=[];\r
3336         for p in cs do\r
3337             mvl:=[op(mvl),Class(p,ord)];\r
3338         od;\r
3339         # notzerolist \r
3340         notzerolist:=CreatNotZeroList(cs,ds,mvl,ord);\r
3341         \r
3342         qv:=[cs,notzerolist[-1]];\r
3343         subasc:=SubASC(cs,pord);\r
3344         lensub:=nops(subasc);\r
3345         rt:=ProjectASC(qv,notzerolist,mvl,subasc,lensub,ord,pord);\r
3346         rt:=Format1(rt);\r
3347 end:\r
3348         \r
3349 \r
3350 #whether input QuasiVariety is [[[],[1]]]\r
3351 QVIsTruth:=proc(qv)\r
3352         if qv=[[[],[1]]] then RETURN(true);\r
3353         else RETURN(false);\r
3354         fi;\r
3355 end:\r
3356 \r
3357 #----main func.-------\r
3358 # project with the QV [ps,ds], eliminate variables in pord\r
3359 # return [] if zero(ps) is empty  \r
3360 # return [[[],[0]]] if zero(ps/ds) is empty\r
3361 # else return [  [[ASC],[p1,p2,...]]  , ...]\r
3362 Project:=proc(ps,ds,ord,pord)\r
3363         local jds,cslist,cs,p,ini,qvs,newds,subasc,rt;\r
3364         jds:=GetProduct(ds);\r
3365         cslist:=wsolve(ps,ord);\r
3366         #printf("The characteristic sets are:");\r
3367         ##lprint(cslist);\r
3368         rt:=[];\r
3369         if nops(cslist)<>0 then\r
3370                 cslist:=Newcslist(cslist,jds,ord);  \r
3371                 if nops(cslist)=0 then RETURN([]);\r
3372                 fi;\r
3373                 cslist:=InverseAll(cslist);\r
3374                 for cs in cslist do\r
3375                         \r
3376                         #add initials to ds\r
3377                         newds:={op(ds)};\r
3378                         for p in cs do\r
3379                                 ini:=Initial(p,ord);\r
3380                                 if Class(ini,ord)>0 then\r
3381                                         newds:={op(newds),ini};\r
3382                                 fi;\r
3383                         od;\r
3384                         newds:=[op(newds)];\r
3385                         \r
3386                         qvs:=ProjASC(cs,newds,ord,pord);\r
3387                         subasc:=qvs[1][1];\r
3388                         qvs:=Simplify1(qvs,subasc,ord);\r
3389                         if QVIsTruth(qvs) then rt:=qvs; RETURN(rt); \r
3390                         else  \r
3391                                 rt:=ListUnion(qvs,rt);\r
3392                         fi;\r
3393                 od;\r
3394         fi;\r
3395         if nops(rt)=0 then RETURN([[[],[0]]]); fi;\r
3396         rt;\r
3397 end:            \r
3398 \r
3399 \r
3400 \r
3401 ## the subasc in qvs is the same one \r
3402 AddQVS:=proc(rt,qvs)\r
3403   local newrt,i,j,subasc;\r
3404   newrt:=rt;\r
3405   subasc:=qvs[1];\r
3406   j:=0;\r
3407   for i from 1 to nops(newrt) do\r
3408      if subasc=newrt[i][1] then newrt[i][2]:=ListUnion(qvs[2],newrt[i][2]); j:=1; break;fi;    \r
3409   od;\r
3410   if j=0 then newrt:=[op(newrt),qvs];fi;\r
3411   newrt;\r
3412 end:\r
3413 \r
3414 SqfreeList:=proc(plist)\r
3415    local j,p,q,rt,sq;\r
3416    rt:=[];\r
3417    for p in plist do\r
3418       j:=1;\r
3419       sq:=sqrfree(p);\r
3420       sq:=sq[2];\r
3421       for q in sq do\r
3422         j:=j*q[1];\r
3423       od;\r
3424       rt:=[op(rt),j];      \r
3425    od;\r
3426    rt;\r
3427 end:\r
3428 \r
3429 \r
3430 SqfreeGB:=proc(plist,ord)\r
3431 local rt,gb1,gb2;\r
3432    gb1:=plist;\r
3433    gb2:=Groebner[gbasis](gb1,plex(op(ord)));\r
3434    gb2:=SqfreeList(gb2);\r
3435    while(gb1<>gb2) do\r
3436       gb1:=gb2;\r
3437       gb2:=Groebner[gbasis](gb1,plex(op(ord)));\r
3438       gb2:=SqfreeList(gb2);   \r
3439    od;\r
3440    gb2;\r
3441 end:    \r
3442 \r
3443 GBsimplify:=proc(qvs,ord)\r
3444 local qv,rt;\r
3445    rt:=[];\r
3446    for qv in qvs do\r
3447        rt:=[op(rt),[qv[1],SqfreeGB(qv[2],ord) ] ];\r
3448    od;\r
3449    rt;\r
3450 end:\r
3451 \r
3452 Congre:=proc(qvs)\r
3453 local rt,qv,notzero;\r
3454     notzero:=[];\r
3455     for qv in qvs do\r
3456         notzero:=[op(notzero),qv[2]];\r
3457      od;\r
3458      notzero:={op(notzero)};\r
3459      notzero:=[op(notzero)];\r
3460     rt:=[qvs[1][1],notzero];\r
3461 end:\r
3462 \r
3463 #new project\r
3464 Proj:=proc(ps,ds,ord,pord)\r
3465         local jds,cslist,cs,p,ini,qvs,newds,subasc,rt;\r
3466         jds:=GetProduct(ds);\r
3467         cslist:=wsolve(ps,ord);\r
3468         #printf("The characteristic sets are:");\r
3469         ##lprint(cslist);\r
3470         rt:=[];\r
3471         if nops(cslist)<>0 then\r
3472                 cslist:=Newcslist(cslist,jds,ord);  \r
3473                 if nops(cslist)=0 then RETURN("Characteristic Set is empty set");\r
3474                 fi;\r
3475                 cslist:=InverseAll(cslist);\r
3476                 for cs in cslist do\r
3477                         \r
3478                         #add initials to ds\r
3479                         newds:={op(ds)};\r
3480                         for p in cs do\r
3481                                 ini:=Initial(p,ord);\r
3482                                 if Class(ini,ord)>0 then\r
3483                                         newds:={op(newds),ini};\r
3484                                 fi;\r
3485                         od;\r
3486                         newds:=[op(newds)];\r
3487                         \r
3488                         qvs:=ProjASC(cs,newds,ord,pord);\r
3489                         subasc:=qvs[1][1];\r
3490                                                 \r
3491                         qvs:=Simplify2(qvs,subasc,ord);\r
3492                         if qvs=[[[],1]] then RETURN([[[],[1]]]); \r
3493                         else  \r
3494                               qvs:=Congre(qvs);\r
3495                               rt:=AddQVS(rt,qvs); \r
3496                         fi;\r
3497                 od;\r
3498         fi;\r
3499         if nops(rt)=0 then RETURN([[[],[0]]]); fi;      \r
3500         rt:=GBsimplify(rt,ord);\r
3501 end:            \r
3502 \r
3503 \r
3504 \r
3505 # return 1 if zero(ps,ds)=empty \r
3506 # else return 0\r
3507 IsEmpty:=proc(ps,ds,ord)\r
3508         local jds,cslist,cs,p,ini,qvs,newds,subasc;\r
3509         if nops(ps)=0 then RETURN(false);fi;\r
3510         jds:=GetProduct(ds);\r
3511         cslist:=wsolve(ps,ord);\r
3512         #printf("The characteristic sets are:");\r
3513         ##lprint(cslist);\r
3514         if nops(cslist)<>0 then\r
3515                 cslist:=Newcslist(cslist,jds,ord);  \r
3516                 if nops(cslist)=0 then RETURN(true);\r
3517                 fi;\r
3518                 cslist:=InverseAll(cslist);\r
3519                 for cs in cslist do\r
3520                         \r
3521                         #add initials to ds\r
3522                         newds:=ds;\r
3523                         for p in cs do\r
3524                                 ini:=Initial(p,ord);\r
3525                                 if Class(ini,ord)>0 then\r
3526                                         newds:=[op(ds),ini];\r
3527                                 fi;\r
3528                         od;\r
3529                         \r
3530                         qvs:=ProjASC(cs,newds,ord,ord);\r
3531                         subasc:=qvs[1][1];\r
3532                         qvs:=Simplify1(qvs,subasc,ord);\r
3533                         if QVIsTruth(qvs) then RETURN(false);\r
3534                         fi;\r
3535                 od;\r
3536         fi;\r
3537         RETURN(true)\r
3538 end:           \r
3539 \r
3540 #return -1 if zero([ ps, [op(ds),p] ])=empty \r
3541 #return 1 if zero([ [op(ps),p],ds ])=empty\r
3542 #else return 0\r
3543 whichcase:=proc(ps0,ds0,p,ord)\r
3544         local ps,ds;\r
3545         ps:=[op(ps0)];\r
3546         ds:=[op(ds0)];\r
3547         if IsEmpty(ps,[op(ds),p],ord,ord) then RETURN(-1) fi;\r
3548         if IsEmpty([op(ps),p],ds,ord,ord) then RETURN(1) fi;            \r
3549         RETURN(0);     \r
3550 end:       \r
3551 \r
3552 \r
3553 #---------------- Format1(qvs)-----------------\r
3554 \r
3555 #when the input is as:["UN",...]\r
3556 #step1: formated all the branchs of "UN".\r
3557 #step2: Union all the branchs.   \r
3558 FormatUN:=proc(qvs)\r
3559         local i,temp,rt;\r
3560         rt:=[];\r
3561         for i from 2 to nops(qvs) do\r
3562                 temp:=Format1(qvs[i]);\r
3563                 rt:=ListUnion(rt,temp);\r
3564         od;\r
3565         rt;\r
3566 end:               \r
3567 \r
3568 # when the input is as:[[ASC],pol]\r
3569 # output [[[ASC],pol]]\r
3570 FormatList:=proc(qv)\r
3571         RETURN([qv]);\r
3572 end:\r
3573 \r
3574 # format the result of ProjectASC()\r
3575 # to the form: [[[ASC],pol1],[[ASC],pol2],...];\r
3576 Format1:=proc(qvs)\r
3577         local rt;\r
3578         if type(qvs[1],string) then\r
3579                 rt:=FormatUN(qvs);\r
3580         else \r
3581                 rt:=FormatList(qvs);\r
3582         fi;\r
3583         rt;\r
3584 end:\r
3585 \r
3586 \r
3587 #-----------------Simplify1(qvs)---------------\r
3588 #for every qausi variety in the result of Project() call doRem1\r
3589 #for empty set , return []; the result form: [[[ASC],[p1,p2,...]],...]\r
3590 #p1<>0 and p2<>0 and ...\r
3591 Simplify1:=proc(qvs,subasc,ord)\r
3592         local qv,rt,newqv;\r
3593         rt:=[];\r
3594         for qv in qvs do\r
3595                 newqv:=doRem1(qv,subasc,ord);\r
3596                 if newqv<>[] then \r
3597                         if member(1,newqv[2]) then RETURN([[subasc,[1]]]);fi;\r
3598                 fi;\r
3599                 rt:=[op(rt),newqv];             \r
3600         od;\r
3601         rt;\r
3602 end:\r
3603 \r
3604 # a qausi variety is denoted as: [[asc],pol]\r
3605 # rem=premas(pol,asc,ord)\r
3606 # notZero:= all factors of rem of dClass w.r.t ord zero.  \r
3607 # return [[asc], [notZero]]\r
3608 doRem1:=proc(qv,subasc,ord)\r
3609         local pol,rem,notZero,rt;\r
3610         pol:=qv[2];\r
3611 # modified   2003.3.21\r
3612 #       if nops(subasc)=0 or Class(pol,ord)=0 then\r
3613 #               if pol=0 then rt:=[];\r
3614 #               else rt:=[subasc,[1]];\r
3615 #               fi;\r
3616 #               RETURN(rt);\r
3617 #       fi;\r
3618         rem:=Premas(pol,Inverse(subasc),ord);\r
3619         if rem=0 then rt:=[];\r
3620         else \r
3621                 notZero:=Factorlist_chen(rem,ord);\r
3622                 rt:=[subasc,notZero];\r
3623         fi;\r
3624         rt;\r
3625 end:\r
3626 #-----------------Simplify2(qvs)---------------\r
3627 #for every qausi variety in the result of Project() call doRem2\r
3628 #for empty set , return []; the result form: [[[ASC],[p1,p2,...]],...]\r
3629 #p1<>0 or p2<>0 or ...\r
3630 \r
3631 #for every qausi variety in the result of Project() call doRem\r
3632 #for empty set , return []; the result form: [[[ASC],[p1,p2,...]],...]\r
3633 Simplify2:=proc(qvs,subasc,ord)\r
3634         local qv,rt,newqv;\r
3635         rt:=[];\r
3636         for qv in qvs do\r
3637                 newqv:=doRem2(qv,subasc,ord);\r
3638                 if newqv<>[] then \r
3639                         if newqv[2]=1 then RETURN([[subasc,1]]);fi;\r
3640                 fi;\r
3641                 rt:=[op(rt),newqv];             \r
3642         od;\r
3643         rt;\r
3644 end:\r
3645 \r
3646 # a qausi variety is denoted as: [[asc],pol]\r
3647 # rem=premas(pol,asc,ord)\r
3648 # return [[asc], rem]\r
3649 doRem2:=proc(qv,subasc,ord)\r
3650         local pol,rem,notZero,rt;\r
3651         pol:=qv[2];\r
3652 \r
3653         rem:=Premas(pol,Inverse(subasc),ord);\r
3654         if rem=0 then rt:=[];\r
3655         else \r
3656                 notZero:=Factorlist_chen(rem,ord);\r
3657                 rt:=[subasc,GetProduct(notZero)];\r
3658         fi;\r
3659         rt;\r
3660 end:\r
3661 \r
3662 \r
3663 # if 0 return 0 else return 1\r
3664 Get1or0:=proc(p)\r
3665         local rt;\r
3666         if p=0 then rt:=0;\r
3667         else rt:=1;\r
3668         fi;\r
3669         rt;\r
3670 end:\r
3671 \r
3672 \r
3673 # p=3*(x+y)^2*(x-y)*u; ord=[x,y]; result=[x+y,x-y];\r
3674 # p=2; result=[1]; p=0; result=[0]\r
3675 Factorlist_chen:=proc(p,ord)\r
3676         local i,j,fl,pair,fact,rt;\r
3677         if Class(p,ord)=0 then rt:=[Get1or0(p)];\r
3678         else \r
3679                 rt:=[];\r
3680                 fl:=factors(p);\r
3681                 fl:=fl[2];\r
3682                 for pair in fl do\r
3683                         fact:=pair[1];\r
3684                         if Class(fact,ord)<>0 then \r
3685                                 rt:=[op(rt),fact];      \r
3686                         fi;\r
3687                 od;\r
3688         fi;\r
3689         rt;\r
3690 end:\r
3691 \r
3692 #########----------Simplify----------###########       \r
3693 \r
3694 \r
3695 ###----- computation of QV---------- ############\r
3696 \r
3697 \r
3698 # judge whether or not Zero(ps/ds) is empty\r
3699 \r
3700 \r
3701 # judge whether or not Zero(ps1/ds1) is included in Zero(ps2/ds2)\r
3702 # QVInclude:=proc(ps1,ds1,ps2,ds2)\r
3703 \r
3704 \r
3705 #------Simplify the output of Project--------# \r
3706 # whether or not QV is included in the union of QVS\r
3707 # output: 1(yes), 0(not included)\r
3708 QVSInclude:=proc(QV,QVS,ord)\r
3709         local D,newQVS,rt,temp,i,ass,proj,ps,ds,newucs;\r
3710         #if QVS=[] then return 0;fi;\r
3711         D:=GetProduct(QV[2]);\r
3712         newQVS:=QVS;\r
3713         for i from 1 to nops(QVS) do \r
3714                 newQVS[i]:=[op(newQVS[i][1]),GetProduct(newQVS[i][2])];         \r
3715         od;\r
3716         ass:=AllSorts(newQVS,ord);\r
3717         for i from 1 to nops(ass) do\r
3718                 temp:=ass[i];\r
3719                 ds:=[D,op(temp[2])];\r
3720                 if temp[1]<>[] then \r
3721                         ps:=[op(QV[1]),op(temp[1])];\r
3722                         proj:=Project(ps,ds,ord,ord);\r
3723                         # case :not empty set \r
3724                         if proj<>[] and proj[1][2][1]<>0 then RETURN(0);fi;\r
3725                 else \r
3726                         proj:=ProjASC(QV[1],ds,ord,ord); \r
3727                         # case: not empty set\r
3728                         if proj[1][2]<>0 then RETURN(0);fi;\r
3729                 fi;\r
3730         od;\r
3731         RETURN(1);\r
3732 end:\r
3733 \r
3734 # if every QV in QVS is included in the union of the others\r
3735 # then expunge it.\r
3736 Sproject:=proc(QVS,ord)\r
3737         local i,rt,UQV,len;\r
3738         rt:=[];\r
3739         len:=nops(QVS);\r
3740         for i from 1 to len do\r
3741                 if i=len then UQV:=[]; fi;\r
3742                 UQV:=QVS[i+1..len];\r
3743                 UQV:=[op(rt),op(UQV)];\r
3744                 if UQV=[] then rt:=[QVS[i]];break ;fi;\r
3745                 if QVSInclude(QVS[i],UQV,ord)=0 then \r
3746                         rt:=[op(rt),QVS[i]];\r
3747                 fi;\r
3748         od;\r
3749         rt;\r
3750 end:\r
3751 \r
3752 Sproj:=proc(qvs,ord)\r
3753         local temprt,rt, tempqvs,qv,nozerop;\r
3754         tempqvs:=[];\r
3755         #change to old form\r
3756         for qv in qvs do \r
3757            for nozerop in qv[2] do\r
3758               tempqvs:=[op(tempqvs), [ qv[1],[nozerop] ] ];\r
3759            od;\r
3760         od;\r
3761         temprt:=Sproject(tempqvs,ord);\r
3762         rt:=[];\r
3763         #return to new form\r
3764         for qv in temprt do\r
3765             rt:=AddQVS(rt,qv);\r
3766         od;\r
3767         rt;\r
3768 end:\r
3769         \r
3770 #------Simplify the output of wsolve--------# \r
3771 #AllSort([[a1,a3],[c1,c2,c3]]);\r
3772 #output:[ [[],[c1,a1]], [[],[c2,a1]], [[c3],[a1]], [[a3],[c1]],\r
3773 #          [[a3],[c2]], [[c3,a3],[]] ]\r
3774 AllSorts:=proc(lists,ord)\r
3775         local ls,C,ls_len,p,i,reculist,list,onesort,rt;\r
3776         \r
3777         rt:=[];\r
3778         ls:=lists[1];\r
3779         ls_len:=nops(lists[1]);\r
3780         \r
3781         C:=Class(ls[ls_len],ord);\r
3782         \r
3783         if nops(lists)=1 then\r
3784                 for i from 1 to ls_len-1 do\r
3785                         p:=[[],[ls[i]]];\r
3786                         rt:=[op(rt),p ];\r
3787                 od;\r
3788                 if C<>0 then\r
3789                         rt:=[op(rt),[[ls[ls_len]],[]]];\r
3790                 fi;\r
3791                 RETURN(rt); \r
3792         fi;\r
3793         \r
3794         reculist:=AllSorts(subsop(1=NULL,lists),ord);\r
3795         \r
3796         for i from 1 to ls_len-1 do\r
3797                 for list in reculist do\r
3798                         onesort:=list;\r
3799                         onesort[2]:=[op(onesort[2]),ls[i]];\r
3800                         rt:=[op(rt),onesort];\r
3801                 od;\r
3802         od;\r
3803         if C<>0 then\r
3804                 for list in reculist do\r
3805                                 onesort:=list;\r
3806                                 onesort[1]:=[op(onesort[1]),ls[ls_len]];\r
3807                                 rt:=[op(rt),onesort];\r
3808                         od;\r
3809         fi;     \r
3810         rt;\r
3811 end:\r
3812 \r
3813 # whether or not cs is included in the union of ucs\r
3814 # output: 1(yes), 0(not included)\r
3815 CSSInclude:=proc(cs,ucs,ord)\r
3816         local proj,IP,ass,rt,temp,i,ps,ds,newucs;\r
3817         if ucs=[] then return 0;fi;\r
3818         IP:=GetIP(cs,ord);\r
3819         newucs:=ucs;\r
3820         for i from 1 to nops(ucs) do \r
3821                 newucs[i]:=[op(newucs[i]),GetIP(newucs[i],ord)];                \r
3822         od;\r
3823         ass:=AllSorts(newucs,ord);\r
3824         for i from 1 to nops(ass) do\r
3825                 temp:=ass[i];\r
3826                 ds:=[IP,op(temp[2])];\r
3827                 if temp[1]<>[] then \r
3828                         ps:=[op(cs),op(temp[1])];\r
3829                         proj:=Project(ps,ds,ord,ord);\r
3830                         # case :not empty set \r
3831                         if proj<>[] and proj[1][2][1]<>0 then RETURN(0);fi;\r
3832                 else \r
3833                         proj:=ProjASC(cs,ds,ord,ord); \r
3834                         # case: not empty set\r
3835                         if proj[1][2]<>0 then RETURN(0);fi;\r
3836                 fi;\r
3837         od;\r
3838         RETURN(1);\r
3839 end:\r
3840 \r
3841 \r
3842 # if every cs in css is included in the union of the others\r
3843 # then expunge it.\r
3844 Swsolve:=proc(css,ord)\r
3845         local i,rt,newcss,UCS,len;\r
3846         rt:=[];\r
3847         newcss:=InverseAll(css);\r
3848         len:=nops(newcss);\r
3849         for i from 1 to len do\r
3850                 UCS:=newcss[i+1..len];\r
3851                 UCS:=[op(rt),op(UCS)];\r
3852                 if CSSInclude(newcss[i],UCS,ord)=0 then \r
3853                         rt:=[op(rt),newcss[i]];\r
3854                 fi;\r
3855         od;\r
3856         InverseAll(rt);\r
3857 end:\r
3858 \r
3859 \r
3860 \r
3861 \r
3862 \r