
# CharSets Version 1.0

######################################################################
#                                                                    #
#                  CHARACTERISTIC SETS PACKAGE                       #
#                                                                    #
#   Author:  Dongming Wang (E-mail: K316170@AEARN.BITNET)            #
#            RISC-LINZ, Johannes Kepler University, Linz, Austria    #
#            Institute of Systems Science, Academia Sinica, China    #
#                                                                    #
#   Date:    December 1990                                           #
#                                                                    #
#   Copyright (C) 1990 by Dongming Wang                              #
#                                                                    #
#   Copyright Notice:  Permission is granted to use, copy or re-     #
#            distribute this package, provided that the title is     #
#            retained and the file is not altered.                   #
#                                                                    #
######################################################################

#====================================================================#
#     This package is implemented for computing characteristic sets  #
#  of (multivariate) polynomial sets,  decomposing  polynomial sets  #
#  into ascending sets and irreducible ascending sets,  decomposing  #
#  algebraic  varieties  into  irreducible components,  factorizing  #
#  polynomials over algebraic number fields  and solving systems of  #
#  polynomial equations.  It is on the basis of  the characteristic  #
#  sets method  introduced by J. F. Ritt  and  developed by Wu Wen-  #
#  tsun. The algorithms with variants implemented here are based on  #
#  a generalization  given by this author.  Other modifications are  #
#  also made. For references, see                                    #
#  Ritt J. F., Differential Algebra, AMS, 1950.                      #
#  Wang D. M., Characteristic Sets and Zero Structure of Polynomial  #
#     Sets, Preprint, RISC-LINZ, 1989.                               #
#  Wu W. T., Basic Principles of Mechanical Theorem Proving in       #
#     Elementary Geometries, J. Sys. Sci. & Math. Scis., 4(1984),    #
#     207-235; J. Automated Reasoning, 2(1986), 221-252.             #
#====================================================================#

##### Part 0. Definition of User Functions #####

charsets[charset] := proc() `charsets/charset`(args) end:

charsets[mcharset] := proc() `charsets/mcharset`(args) end:

charsets[charser] := proc() `charsets/charser`(args) end:

charsets[mcs] := proc() `charsets/mcs`(args) end:

charsets[ecs] := proc() `charsets/ecs`(args) end:

charsets[mecs] := proc() `charsets/mecs`(args) end:

charsets[ics] := proc() `charsets/ics`(args) end:

charsets[qics] := proc() `charsets/qics`(args) end:

charsets[eics] := proc() `charsets/eics`(args) end:

charsets[ivd] := proc() `charsets/ivd`(args) end:

charsets[remset] := proc() `charsets/remset`(args) end:

charsets[cfactor] := proc() `charsets/cfactor`(args) end:

charsets[iniset] := proc() `charsets/iniset`(args) end:

charsets[csolve] := proc() `charsets/csolve`(args) end:

charsets[triser] := proc() `charsets/triser`(args) end:

##### Part I. Routines for Computing Characteristic Sets #####

# the class of poly f wrt variable ordering ord
`charsets/class` := proc(f,ord)
                  local i;
                  options remember,system;
                      for i from nops(ord) by -1 to 1 do
                          if has(f,ord[i]) then RETURN(i) fi
                      od;
                      0
                  end:

# the leading variable of poly f wrt variable ordering ord
`charsets/lvar` := proc(f,ord) ord[`charsets/class`(f,ord)] end:

# the index set of a poly (or a poly set f) wrt ord
`charsets/index` :=

    proc(f,ord)
    local i;
        if type(f,list) then
            [seq(`charsets/index`(f[i],ord), i = 1 .. nops(f))]
        elif type(f,set) then
            {seq(`charsets/index`(f[i],ord), i = 1 .. nops(f))}
        else
            if `charsets/class`(f,ord) = 0 then [nops(expand(f)),0,0]
            else
                [nops(expand(f)),`charsets/class`(f,ord),
                    degree(f,`charsets/lvar`(f,ord))]
            fi
        fi
    end:

# the initial of poly f wrt ord
`charsets/initial` :=

    proc(f,ord)
    options remember,system;
        if `charsets/class`(f,ord) = 0 then 1
        else lcoeff(f,`charsets/lvar`(f,ord)); numer("/lcoeff("))
        fi
    end:

# modified rank of two polys: comparing further the rank
#     of initials when f and g have same rank
`charsets/mrank` :=

proc(f,g,ord)
local cf,cg;
options remember,system;
    cf := `charsets/class`(f,ord);
    cg := `charsets/class`(g,ord);
    if cf = 0 then true
    elif cf < cg then true
    elif cf = cg then
        cf := degree(f,`charsets/lvar`(f,ord));
        cg := degree(g,`charsets/lvar`(g,ord));
        if cf < cg then true
        elif cf = cg then
            `charsets/mrank`(
                `charsets/initial`(f,ord),`charsets/initial`(g,ord),ord)
        else false
        fi
    else false
    fi
end:

# modified rank of two polys: comparing further the rank of
#     initials, the terms of initials and the terms of f and g
#     when they are the same
`charsets/rank` :=

    proc(f,g,ord)
    local ind,find,cf,cg;
    options remember,system;
        find := `charsets/subrank`(f,g,ord,'ind');
        if find and ind = 1 then
            cf := nops(expand(`charsets/initial`(f,ord)));
            cg := nops(expand(`charsets/initial`(g,ord)));
            if cf < cg then true
            elif cf = cg then
                if nops(expand(f)) < nops(expand(g)) then true else false fi
            else false
            fi
        else find
        fi
    end:

# subroutine for rank
`charsets/subrank` :=

proc(f,g,ord,ind)
local cf,cg;
options remember,system;
    cf := `charsets/class`(f,ord);
    cg := `charsets/class`(g,ord);
    if cf = 0 then
        if cg = 0 then ind := 1 fi; true
    elif cf < cg then true
    elif cf = cg then
        cf := degree(f,`charsets/lvar`(f,ord));
        cg := degree(g,`charsets/lvar`(g,ord));
        if cf < cg then true
        elif cf = cg then
            `charsets/subrank`(
                `charsets/initial`(f,ord),`charsets/initial`(g,ord),ord,'ind')
        else false
        fi
    else false
    fi
end:

# the rank of two polys with same classes:
#        used for computing tiangular form
`charsets/trank` :=

proc(f,g,ord)
local cf,cg;
options remember,system;
    cf := degree(f,`charsets/lvar`(f,ord));
    cg := degree(g,`charsets/lvar`(g,ord));
    if cf < cg then true
    elif cf = cg then
        `charsets/mrank`(
            `charsets/initial`(f,ord),`charsets/initial`(g,ord),ord)
    else false
    fi
end:

# modified pseudo division: I1^s1...Ir^sr*uu = q*vv + r,
#    where I1, ..., I_r are all distinct factors of lcoeff(vv,x)
#    and s1, ..., sr are chosen to be the smallest
`charsets/prem` :=

    proc(uu,vv,x)
    local r,v,dr,dv,l,t,lu,lv;
    options remember,system;
        if type(vv/x,integer) then subs(x = 0,uu)
        else
            r := expand(uu);
            dr := degree(r,x);
            v := expand(vv);
            dv := degree(v,x);
            if dv <= dr then l := coeff(v,x,dv); v := expand(v-l*x^dv)
            else l := 1
            fi;
            while dv <= dr and r <> 0 do
                gcd(l,coeff(r,x,dr),'lu','lv');
                t := expand(x^(dr-dv)*v*lv);
                if dr = 0 then r := 0 else r := subs(x^dr = 0,r) fi;
                r := expand(lu*r)-t;
                dr := degree(r,x)
            od;
            r
        fi
    end:

# pseudo remainder of poly f wrt ascending set as
`charsets/premas` :=

    proc(f,as,ord)
    local remd,i;
        remd := f;
        for i from nops(as) by -1 to 1 do
            remd := `charsets/prem`(remd,as[i],`charsets/lvar`(as[i],ord))
        od;
        if remd <> 0 then numer(remd/lcoeff(remd)) else 0 fi
    end:

# set of non-zero remainders of polys in ps wrt ascending set as
#       user level function
`charsets/remset` :=

proc(ps,as,ord)
local ind,i;
    if nargs <> 3 then ERROR(`wrong number of arguments`)
    elif nops(ps) < 1 or nops(as) < 1 then ERROR(`no polynomials specified`)
    elif nops(ord) < 1 then ERROR(`no indeterminates specified`)
    elif not type(ord,list) then ERROR(ord,`must be a list`)
    fi;
    if member(false,map(type,ord,name)) then ERROR(`bad variable list`) fi;
    ind := 0;
    for i to nops(as) do
        if `charsets/class`(as[i],ord) <= ind then
            ERROR(
    `second argument must be a non-contradictory (weak, quasi-) ascending set`
            )
        else ind := `charsets/class`(as[i],ord)
        fi
    od;
    if type(ps,{set,list}) then
        if member(false,map(type,ps,polynom(rational,ord))) or
            member(false,map(type,as,polynom(rational,ord))) then
            ERROR(`input must be polynomials over Q in`,ord)
        fi;
        `charsets/remseta`(ps,as,ord)
    else
        if member(false,map(type,as,polynom(rational,ord))) then
            ERROR(`input must be polynomials over Q in`,ord)
        fi;
        `charsets/premas`(ps,as,ord)
    fi
end:

# set of non-zero remainders of polys in ps wrt ascending set as
`charsets/remseta` :=

    proc(ps,as,ord)
    local i;
        {seq(`charsets/premas`(ps[i],as,ord), i = 1 .. nops(ps))} minus {0}
    end:

# pseudo remainder of poly f wrt ascending set as -- version b
`charsets/premasb` :=

    proc(f,as,ord)
    local remd,i;
        remd := f;
        if nops(as) > 1 then
            for i from nops(as) by -1 to 2 do
                remd := `charsets/prem`(remd,as[i],`charsets/lvar`(as[i],ord))
            od
        fi;
        if divide(remd,as[1]) then remd := 0
        else remd := `charsets/prem`(remd,as[1],`charsets/lvar`(as[1],ord)) fi;
        if remd <> 0 then numer(remd/lcoeff(remd)) else 0 fi
    end:

# set of non-zero remainders of polys in ps wrt ascending set as -- version b
`charsets/remsetb` :=

    proc(ps,as,ord)
    local i;
        {seq(`charsets/premasb`(ps[i],as,ord), i = 1 .. nops(ps))} minus {0}
    end:

# reorder the list ord of variables wrt polyset ps
`charsets/reorder` :=

    proc(ord,p,ps)
        op(`charsets/reordera`(ord,ps));
        [op(`charsets/reorderb`([op({op(ord)} minus {"})],p,ps)),"]
    end:

# subroutine for reorder: first criterion
`charsets/reordera` :=

proc(ord,ps)
local qs,pp,orb,i;
    if nops(ps) = 0 then ord
    else
        qs := {op(ps)};
        orb := {op(ord)};
        for i in orb do
            pp := `charsets/deg0`(ps,i);
            if nops(pp) = 1 then
                RETURN(
                 [op(`charsets/reordera`([op(orb minus {i})],qs minus pp)),i]
                 )
            fi
        od;
        []
    fi
end:

# subroutine for reorder -- modified from sort: second criterion
`charsets/reorderb` :=

    proc(l,p,ps)
    local n,tn,gap,i,j,temp,v;
        n := nops(l);
        tn := p;
        for i to n do  v[i-1] := l[i] od;
        for gap from 4 while gap <= n do  gap := 3*gap+1 od;
        gap := iquo(gap,3);
        while 0 < gap do
            for i from gap to n-1 do
                temp := v[i];
                for j from i-gap by -gap to 0 do
                    if tn(v[j],temp,ps) then break fi; v[j+gap] := v[j]
                od;
                v[j+gap] := temp
            od;
            gap := iquo(gap,3)
        od;
        [seq( v[i], i=0..n-1 )]
    end:

# determine the order between x and y wrt ps
`charsets/degord` :=

proc(x,y,ps)
    if op(2,`charsets/degpsmax`(ps,y)) < op(2,`charsets/degpsmax`(ps,x)) then
        true
    elif op(2,`charsets/degpsmax`(ps,x)) < op(2,`charsets/degpsmax`(ps,y)) then
        false
    elif op(1,`charsets/degpsmax`(ps,y)) < op(1,`charsets/degpsmax`(ps,x)) then
        true
    elif op(1,`charsets/degpsmax`(ps,x)) < op(1,`charsets/degpsmax`(ps,y)) then
        false
    elif op(2,`charsets/degpsmin`(ps,x)) < op(2,`charsets/degpsmin`(ps,y)) then
        true
    elif op(2,`charsets/degpsmin`(ps,y)) < op(2,`charsets/degpsmin`(ps,x)) then
        false
    elif op(1,`charsets/degpsmin`(ps,y)) < op(1,`charsets/degpsmin`(ps,x)) then
        true
    elif op(1,`charsets/degpsmin`(ps,x)) < op(1,`charsets/degpsmin`(ps,y)) then
        false
    elif op(1,`charsets/deg1`(ps,y)) < op(1,`charsets/deg1`(ps,x)) then true
    elif op(1,`charsets/deg1`(ps,x)) < op(1,`charsets/deg1`(ps,y)) then false
    elif op(2,`charsets/deg1`(ps,y)) < op(2,`charsets/deg1`(ps,x)) then true
    else false
    fi
end:

# the maximal degree of polys in ps wrt x
#      and the number of polys having this degree
`charsets/degpsmax` :=

    proc(ps,x)
    local i,m,mm;
    options remember,system;
        m := max(seq(degree(ps[i],x), i = 1 .. nops(ps)));
        mm := 0;
        for i to nops(ps) do  if degree(ps[i],x) = m then mm := mm+m fi od;
        [mm,m]
    end:

# the minimal non-zero degree of polys in ps wrt x
#      and the number of polys having this degree
`charsets/degpsmin` :=

    proc(ps,x)
    local i,m,mm;
    options remember,system;
        {seq(degree(ps[i],x), i = 1 .. nops(ps))} minus {0};
        if " = {} then m := 0 else m := min(op(")) fi;
        mm := 0;
        for i to nops(ps) do  if degree(ps[i],x) = m then mm := mm+m fi od;
        [mm,m]
    end:

# determine if ps has one and only one poly involving x
`charsets/deg0` := proc(ps,x)
                 local i,ms;
                     ms := {};
                     for i in ps while nops(ms) < 2 do
                         if has(i,x) then ms := {op(ms),i} fi
                     od;
                     ms
                 end:

# the minimal total degree of lcoeffs of polys in ps wrt x
#      and the minimal number of terms of those lcoeffs
`charsets/deg1` :=

    proc(ps,x)
    local i,qs,k;
    options remember,system;
        qs := {};
        k := op(2,`charsets/degpsmin`(ps,x));
        for i to nops(ps) do
            if degree(ps[i],x) = k then qs := {op(qs),lcoeff(ps[i],x)} fi
        od;
        [min(seq(degree(qs[i],indets(qs[i])), i = 1 .. nops(qs))),
            min(seq(nops(expand(qs[i])), i = 1 .. nops(qs)))]
    end:

# search an element with lowest rank in ps
#      and assign the rest of polys to qs
`charsets/sort` :=

    proc(ps,rank,ord,qs)
    local l,i,qs1;
        if nops(ps) = 1 then qs := []; ps[1]
        else
            l := ps[1];
            qs1 := [];
            for i from 2 to nops(ps) do
                if rank(ps[i],l,ord) then qs1 := [l,op(qs1)]; l := ps[i]
                else qs1 := [ps[i],op(qs1)]
                fi
            od;
            qs := qs1;
            l
        fi
    end:

# the difference of two lists
`charsets/minus` := proc(ps,qs) [op({op(ps)} minus {op(qs)})] end:

# the union of three lists
`charsets/union` :=

    proc(ps1,ps2,ps3) [op(({op(ps1)} union {op(ps2)}) union {op(ps3)})] end:

# the product of all elements in a list
`charsets/prod` := proc(ps) local i; product('ps[i]','i' = 1..nops(ps)) end:

# the basic set of polyset ps
`charsets/basset` := proc(ps,ord)
                   local qs,qs1,i,b;
                       if nops(ps) < 2 then ps
                       else
                           b := `charsets/sort`(ps,`charsets/rank`,ord,'qs1');
                           qs := [];
                           if 0 < `charsets/class`(b,ord) then
                               for i in qs1 do
                                   if degree(i,`charsets/lvar`(b,ord)) <
                                       degree(b,`charsets/lvar`(b,ord)) then
                                       qs := [i,op(qs)]
                                   fi
                               od
                           else RETURN([b])
                           fi;
                           [b,op(`charsets/basset`(qs,ord))]
                       fi
                   end:

# the weak basic set of polyset ps
`charsets/wbasset` :=

proc(ps,ord)
local qs,qs1,i,b;
    if nops(ps) < 2 then ps
    else
        b := `charsets/sort`(ps,`charsets/rank`,ord,'qs1');
        qs := [];
        if 0 < `charsets/class`(b,ord) then
            for i in qs1 do
                if `charsets/class`(b,ord) < `charsets/class`(i,ord) and
                    degree(`charsets/initial`(i,ord),`charsets/lvar`(b,ord)) <
                    degree(b,`charsets/lvar`(b,ord)) then
                    qs := [i,op(qs)]
                fi
            od
        else RETURN([b])
        fi;
        [b,op(`charsets/wbasset`(qs,ord))]
    fi
end:

# the quasi-basic set of polyset ps
`charsets/qbasset` :=

    proc(ps,ord)
    local qs,qs1,i,b;
        if nops(ps) < 2 then ps
        else
            b := `charsets/sort`(ps,`charsets/rank`,ord,'qs1');
            qs := [];
            if 0 < `charsets/class`(b,ord) then
                for i in qs1 do
                    if `charsets/class`(b,ord) < `charsets/class`(i,ord) then
                        qs := [i,op(qs)]
                    fi
                od
            else RETURN([b])
            fi;
            [b,op(`charsets/qbasset`(qs,ord))]
        fi
    end:

# the char set of polyset ps: user function
`charsets/charset` :=

proc(ps,lst,medset,y)
local mset,ord;
    if nargs < 2 then ERROR(`too few arguments`)
    elif nops(ps) < 1 then ERROR(`no polynomials specified`)
    elif nops(lst) < 1 then ERROR(`no indeterminates specified`)
    elif 4 < nargs then ERROR(`too many arguments`)
    fi;
    if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi;
    if member(false,map(type,[op(ps)],polynom(rational,lst))) then
        ERROR(`input must be polynomials over Q in`,lst)
    fi;
    if type(lst,list) then ord := lst
    else ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps)
    fi;
    if nargs < 3 then mset := 'charsetn' else mset := medset fi;
    if 3 < nargs then y := ord fi;
    if member(
      mset,{'wcharsetn','charsetn','qcharsetn','wbasset','qbasset','triset',
            'trisetc','basset'}) then
        `charsets/charseta`(expand(ps),ord,`charsets/`.mset)
    else
        ERROR(`medial set must be one of ``basset``,``wbasset``, ``qbasset``,`.
        ```charsetn``,``wcharsetn``,``qcharsetn``,``triset`` and ``trisetc```)
    fi
end:

# the char set of polyset ps
`charsets/charseta` :=

    proc(ps,ord,medset)
    local cs,rs,l;
        if nops(ps) < 2 then ps
        else
            cs := medset(ps,ord);
            if 0 < `charsets/class`(cs[1],ord) then
                if member(medset,{`charsets/basset`,`charsets/wbasset`,
                                  `charsets/qbasset`}) then
                    rs := `charsets/remseta`({op(ps)} minus {op(cs)},cs,ord)
                else rs := `charsets/remsetb`({op(ps)} minus {op(cs)},cs,ord)
                fi
            else RETURN([1])
            fi;
            if rs = {} then
                [seq(numer(cs[l]/lcoeff(cs[l])), l = 1 .. nops(cs))]
            else `charsets/charseta`(`charsets/union`(rs,cs,ps),ord,medset)
            fi
        fi
    end:

# the modified char set of polyset ps: user function
`charsets/mcharset` :=

proc(ps,lst,medset,y)
local mset,ord;
    if nargs < 2 then ERROR(`too few arguments`)
    elif nops(ps) < 1 then ERROR(`no polynomials specified`)
    elif nops(lst) < 1 then ERROR(`no indeterminates specified`)
    elif 4 < nargs then ERROR(`too many arguments`)
    fi;
    if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi;
    if member(false,map(type,[op(ps)],polynom(rational,lst))) then
        ERROR(`input must be polynomials over Q in`,lst)
    fi;
    if type(lst,list) then ord := lst
    else ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps)
    fi;
    if 3 < nargs then y := ord fi;
    if nargs < 3 then mset := 'charsetn' else mset := medset fi;
    if member(
      mset,{'wcharsetn','charsetn','qcharsetn','wbasset','qbasset','triset',
            'trisetc','basset'}) then
        `charsets/fcharseta`(expand(ps),ord,`charsets/`.mset)
    else
        ERROR(`medial set must be one of ``basset``,``wbasset``, ``qbasset``,`.
        ```charsetn``,``wcharsetn``,``qcharsetn``,``triset`` and ``trisetc```)
    fi
end:

# the modified char set of polyset ps
`charsets/fcharseta` :=

    proc(ps,ord,medset)
    local csf,i,fset;
        csf := `charsets/fcharsetsub`(`charsets/nopower`(ps),ord,medset,
            [{},{seq(op(indets(ps[i])), i = 1 .. nops(ps))}],'fset');
        csf,`factors removed` = fset[1]
    end:

# the main subroutine for fcharseta
`charsets/fcharsetsub` :=

 proc(ps,ord,medset,fset1,fset)
 local cs,rs,l,fset3,fset2,ts,fmedset;
     if nops(ps) < 2 then fset := fset1; ps
     else
         if member(substring(medset,10 .. length(medset)),
             {'wcharsetn','charsetn','qcharsetn','triset','trisetc'}) then
             fmedset := `charsets/f`.(substring(medset,10 .. length(medset)));
             cs := fmedset(ps,ord,fset1,'fset2');
             if nops(indets(cs[1])) > 2 then
                 cs := `charsets/removecont`(cs,ord,'ts');
                 fset2 := [fset2[1] union ts,fset2[2]]
             fi;
             if 0 < `charsets/class`(cs[1],ord) then
                 rs := `charsets/remsetb`({op(ps)} minus {op(cs)},cs,ord)
             else fset := fset2; RETURN([1])
             fi;
             if rs = {} then
                 fset := fset2;
                 [seq(numer(cs[l]/lcoeff(cs[l])), l = 1 .. nops(cs))]
             else
                 `charsets/fcharsetsub`(
                     `charsets/union`(rs,cs,ps),ord,medset,fset2,'fset')
             fi
         else
             cs := medset(ps,ord);
             fset2 := [fset1[1],fset1[2] union `charsets/initialset1`(cs,ord)];
             if 0 < `charsets/class`(cs[1],ord) then
                 `charsets/remseta`({op(ps)} minus {op(cs)},cs,ord);
                 rs := `charsets/removefactor`(",ord,fset2,'fset3')
             else fset := fset2; RETURN([1])
             fi;
             if rs = [] then
                 fset := fset3;
                 [seq(numer(cs[l]/lcoeff(cs[l])), l = 1 .. nops(cs))]
             else
                 `charsets/fcharsetsub`(
                     `charsets/union`(rs,cs,ps),ord,medset,fset3,'fset')
             fi
         fi
     fi
 end:

# replace the power of factors of polys in as by 1 if any
`charsets/nopower` :=

  proc(as)
  local i;
      if not type(as,{set,list}) then
          if type(as,`^`) then op(1,as)
          elif type(as,`*`) then
              product('`charsets/nopower`(op(i,as))','i' = 1 .. nops(as))
          else as
          fi
      else
          [seq(`charsets/nopower`(as[i]), i = 1 .. nops(as))]
      fi
  end:

# the nearly char set -- a medial set
`charsets/charsetn` :=

    proc(ps,ord)
    local cs,rs;
        if nops(ps) < 2 then ps
        else
            cs := `charsets/basset`(ps,ord);
            if 0 < `charsets/class`(cs[1],ord) then
                rs := `charsets/remseta`(`charsets/minus`(ps,cs),cs,ord)
            else RETURN(cs)
            fi;
            if rs = {} then cs else `charsets/charsetn`([op(rs),op(cs)],ord) fi
        fi
    end:

# the nearly weak char set -- a weak medial set
`charsets/wcharsetn` :=

   proc(ps,ord)
   local cs,rs;
       if nops(ps) < 2 then ps
       else
           cs := `charsets/wbasset`(ps,ord);
           if 0 < `charsets/class`(cs[1],ord) then
               rs := `charsets/remseta`(`charsets/minus`(ps,cs),cs,ord)
           else RETURN(cs)
           fi;
           if rs = {} then cs else `charsets/wcharsetn`([op(rs),op(cs)],ord) fi
       fi
   end:

# the nearly quasi-char set -- a quasi-medial set
`charsets/qcharsetn` :=

   proc(ps,ord)
   local cs,rs;
       if nops(ps) < 2 then ps
       else
           cs := `charsets/qbasset`(ps,ord);
           if 0 < `charsets/class`(cs[1],ord) then
               rs := `charsets/remseta`(`charsets/minus`(ps,cs),cs,ord)
           else RETURN(cs)
           fi;
           if rs = {} then cs else `charsets/qcharsetn`([op(rs),op(cs)],ord) fi
       fi
   end:

# the modified nearly char set -- a modified medial set
`charsets/fcharsetn` :=

    proc(ps,ord,fset1,fset)
    local cs,rs,fset2,fset3;
        if nops(ps) < 2 then fset := fset1; ps
        else
            cs := `charsets/basset`(ps,ord);
            fset2 := [fset1[1],fset1[2] union `charsets/initialset1`(cs,ord)];
            if 0 < `charsets/class`(cs[1],ord) then
                `charsets/remseta`(`charsets/minus`(ps,cs),cs,ord);
                rs := `charsets/removefactor`(",ord,fset2,'fset3')
            else fset := fset2; RETURN(cs)
            fi;
            if rs = [] then fset := fset3; cs
            else `charsets/fcharsetn`([op(rs),op(cs)],ord,fset3,'fset')
            fi
        fi
    end:

# the modified nearly weak char set -- a modified weak medial set
`charsets/fwcharsetn` :=

    proc(ps,ord,fset1,fset)
    local cs,rs,fset2,fset3;
        if nops(ps) < 2 then fset := fset1; ps
        else
            cs := `charsets/wbasset`(ps,ord);
            fset2 := [fset1[1],fset1[2] union `charsets/initialset1`(cs,ord)];
            if 0 < `charsets/class`(cs[1],ord) then
                `charsets/remseta`(`charsets/minus`(ps,cs),cs,ord);
                rs := `charsets/removefactor`(",ord,fset2,'fset3')
            else fset := fset2; RETURN(cs)
            fi;
            if rs = [] then fset := fset3; cs
            else `charsets/fwcharsetn`([op(rs),op(cs)],ord,fset3,'fset')
            fi
        fi
    end:

# the modified nearly quasi-char set -- a modified quasi-medial set
`charsets/fqcharsetn` :=

    proc(ps,ord,fset1,fset)
    local cs,rs,fset2,fset3;
        if nops(ps) < 2 then fset := fset1; ps
        else
            cs := `charsets/qbasset`(ps,ord);
            fset2 := [fset1[1],fset1[2] union `charsets/initialset1`(cs,ord)];
            if 0 < `charsets/class`(cs[1],ord) then
                `charsets/remseta`(`charsets/minus`(ps,cs),cs,ord);
                rs := `charsets/removefactor`(",ord,fset2,'fset3')
            else fset := fset2; RETURN(cs)
            fi;
            if rs = [] then fset := fset3; cs
            else `charsets/fqcharsetn`([op(rs),op(cs)],ord,fset3,'fset')
            fi
        fi
    end:

# the triangular set of polyset ps -- a quasi-medial set
`charsets/triset` :=

   proc(ps,ord)
   local i;
       if nops(ps) < 2 then ps
       else
           for i from 0 to nops(ord) do  `charsets/@qs`[i] := [] od;
           for i in ps do
               `charsets/@qs`[`charsets/class`(i,ord)] :=
                   [op(`charsets/@qs`[`charsets/class`(i,ord)]),i]
           od;
           for i from nops(ord) by -1 to 1 do
               if `charsets/@qs`[0] = [] then
                   `charsets/subtriset`(i,ord)
               else RETURN([1])
               fi
           od;
           if `charsets/@qs`[0] <> [] then [1]
           else
               [seq(op(`charsets/@qs`[i]), i = 1 .. nops(ord))]
           fi
       fi
   end:

# subroutine for triset
`charsets/subtriset` :=

proc(i,ord)
local ss,ss1,j,p;
    if 1 < nops(`charsets/@qs`[i]) then
        ss1 := `charsets/sort`(`charsets/@qs`[i],`charsets/trank`,
               ord,'ss');
        `charsets/@qs`[i] := [ss1];
        for j in ss do
            p := `charsets/prem`(j,ss1,ord[i]);
            if p <> 0 then
                `charsets/@qs`[`charsets/class`(p,ord)] :=
                 [op(`charsets/@qs`[`charsets/class`(p,ord)]),
                  numer(p/lcoeff(p))]
            fi
        od;
        `charsets/subtriset`(i,ord)
    fi
end:

# remove factor g from f until f has no factor g
#      if g is removed then assign true to ja
`charsets/movefactor` :=

    proc(f,g,ord,ja)
    local fg;
        if (not type(g,integer)) and divide(f,g,'fg') then
            if 3 < nargs then
                if 0 < `charsets/class`(g,ord) then ja := true else
                    ja := false fi
            fi;
            `charsets/movefactor`(fg,g,ord)
        else if 3 < nargs then ja := false fi; f
        fi
    end:

# remove possible factors in fset1[1] and fset1[2] from polys in ps
#      where fset1[1] contains all factors removed before
#      if any poly in fset1[2] is removed, it is added to fset1[1]
#      fset1 is assigned to fset at th end of the procedure
`charsets/removefactor` :=

    proc(ps,ord,fset1,fset)
    local k,rr,ja,fset2,fs,qs,rs,r;
        if not type(ps,{set,list}) then qs := {ps} else qs := ps fi;
        rs := {};
        fs := fset1;
        for r in qs do
            rr := r;
            fset2 := {};
            for k in fs[1] do
                rr := `charsets/movefactor`(rr,k,ord,'ja')
            od;
            for k in fs[2] do
                if rr <> k then
                    rr := `charsets/movefactor`(rr,k,ord,'ja');
                    if ja then fset2 := {numer(k/lcoeff(k)),op(fset2)} fi
                fi
            od;
            fs := [fs[1] union fset2,fs[2] minus fset2];
            rs := {rr,op(rs)}
        od;
        rs := `charsets/nopower`(rs);
        fset := fs;
        if not type(ps,{set,list}) then rs[1] else rs fi
    end:

# remove contents of all polys in ps wrt leading variables
#      the set of removed factors is assigned to ms
`charsets/removecont` :=

  proc(ps,ord,ms)
  local qs,fs,i,cc,pp;
      if `charsets/class`(ps[1],ord) = 0 then if 2 < nargs then ms := {} fi; ps
      else
          qs := [];
          fs := {};
          for i to nops(ps) do
              cc := content(ps[i],`charsets/lvar`(ps[i],ord),'pp');
              if 0 < `charsets/class`(cc,ord) then fs := {op(fs),cc} fi;
              qs := [op(qs),pp]
          od;
          if 2 < nargs then ms := fs fi;
          qs
      fi
  end:

# the modified triangular set of polyset ps -- a modified quasi-medial set
`charsets/ftriset` :=

proc(ps,ord,fset1,fset)
local i,fset2,var;
    `charsets/@fact` := 1;
    if nops(ps) < 2 then fset := fset1; ps
    else
        fset2 := {};
        var := {seq(op(indets(ps[i])), i = 1 .. nops(ps))};
        for i from 0 to nops(ord) do  `charsets/@qs`[i] := [] od;
        for i in ps do
            `charsets/@qs`[`charsets/class`(i,ord)] :=
                [op(`charsets/@qs`[`charsets/class`(i,ord)]),i]
        od;
        for i from nops(ord) by -1 to 1 do
            if `charsets/@qs`[0] = [] then
                fset2 :=
                 fset2 union `charsets/fsubtriset`(i,ord,ps,var)
            else fset := [fset2,{}]; RETURN([1])
            fi
        od;
        if `charsets/@qs`[0] <> [] then fset := [fset2,{}]; [1]
        else
            fset := [fset2,{}];
            [seq(op(`charsets/@qs`[i]), i = 1 .. nops(ord))]
        fi
    fi
end:

# subroutine for ftriset
`charsets/fsubtriset` :=

  proc(i,ord,ps,var)
  local fset2,ss,ss1,j,k,l,p,pp,qq,ja;
      if 1 < nops(`charsets/@qs`[i]) then
          ss1 := `charsets/sort`(`charsets/@qs`[i],`charsets/trank`,ord,'ss');
          `charsets/@qs`[i] := [ss1];
          fset2 := {};
          qq := {seq(numer(ps[l]/lcoeff(ps[l])), l = 1 .. nops(ps))};
          for j in ss do
              pp := `charsets/prem`(j,ss1,ord[i]);
              if pp <> 0 then
                  if pp <> `charsets/@fact` and
                      `charsets/fsubtrisetsub`(var,`charsets/@fact`) and
                      (not member(`charsets/@fact`,qq)) then
                      p := `charsets/movefactor`(pp,`charsets/@fact`,ord,'ja');
                      if ja then
                          fset2 := {op(fset2),
                              numer(`charsets/@fact`/lcoeff(`charsets/@fact`))}
                      fi
                  else p := pp
                  fi;
                  for k in var do
                      if p <> k and (not member(k,qq)) then
                          p := `charsets/movefactor`(p,k,ord,'ja');
                          if ja then fset2 := {op(fset2),k} fi
                      fi
                  od;
                  `charsets/@qs`[`charsets/class`(p,ord)] := [
                      op(`charsets/@qs`[`charsets/class`(p,ord)]),
                      `charsets/nopower`(numer(p/lcoeff(p)))]
              fi
          od;
          `charsets/initial`(ss1,ord);
          `charsets/@fact` := numer("/lcoeff("));
          fset2 union `charsets/fsubtriset`(i,ord,ps,var)
      else {}
      fi
  end:

# subroutine for fsubtriset
`charsets/fsubtrisetsub` :=

  proc(aa,bb)
  local i;
      for i to nops(aa) do  if subs(aa[i] = 0,bb) = 0 then RETURN(false) fi od;
      true
  end:

# reduce a triangular set into an ascending set
`charsets/trisetc` :=

    proc(ps,ord)
    local ind,cs;
        `charsets/@cs` := `charsets/triset`(ps,ord);
        cs := `charsets/subtrisetc`(ord,{},'ind');
        if ind = 0 then `charsets/charsetn`([op(ps),op(cs)],ord) else cs fi
    end:

# reduce a triangular set into an ascending set with factors moved
`charsets/ftrisetc` :=

    proc(ps,ord,fset1,fset)
    local i,ind,cs,fs;
        `charsets/@cs` := `charsets/ftriset`(ps,ord,fset1,'fs');
        fset := fs;
        if fs[1] <> {} then
            {seq(op(`charsets/pfactor`(fs[1][i])), i = 1 .. nops(fs[1]))}
        else {}
        fi;
        cs := `charsets/subtrisetc`(ord,",'ind');
        if ind = 0 then `charsets/fcharsetn`([op(ps),op(cs)],ord,fset1,'fset')
        else cs
        fi
    end:

# subroutine for ftrisetc
`charsets/subtrisetc` :=

proc(ord,var,ind)
local r,i,j,cs;
    if nops(`charsets/@cs`) = 0 then cs := []
    else cs := [`charsets/@cs`[1]]
    fi;
    if 1 < nops(`charsets/@cs`) then
        for i from 2 to nops(`charsets/@cs`) do
            r := `charsets/premas`(`charsets/@cs`[i],cs,ord);
            if
             `charsets/class`(r,ord) <> `charsets/class`(`charsets/@cs`[i],ord)
              then
                ind := 0;
                if r <> 0 then
                    cs := [op(cs),`charsets/nopower`(numer(r/lcoeff(r)))]
                fi;
                break
            else
                for j in var do  r := `charsets/movefactor`(r,j,ord) od;
                cs := [op(cs),`charsets/nopower`(numer(r/lcoeff(r)))];
                if `charsets/class`(r,ord) <>
                    `charsets/class`(`charsets/@cs`[i],ord) then
                    ind := 0; break
                fi
            fi
        od
    fi;
    cs
end:

# the set of nonconstant initials of as
#    with certain repeated factors cancelled
`charsets/initialset1` :=

proc(as,ord)
    local i,is,iss;
    is := {};
    for i in as do
        `charsets/initial`(i,ord);
        if `charsets/class`(",ord) > 0 then
            is := {op(is),`charsets/pfactor`(")}
        fi
    od;
    is := `charsets/compress`(is,ord);
    iss := {};
    for i in is do
        if `charsets/class`(i,ord) > 0 then iss := {op(iss),i} fi
    od;
    iss
end:

# compress some repeated factors
`charsets/compress` :=

    proc(ps,ord)
    local is,is1,i,j,ss;
        is := ps;
        if 1 < nops(is) then
            is1 := [];
            for i to nops(is)-1 do
                ss := is[i];
                for j from i+1 to nops(is) do
                    ss := `charsets/movefactor`(ss,is[j],ord)
                od;
                is1 := [`charsets/pfactor`(ss),op(is1)]
            od;
            is1 := [is[nops(is)],op(is1)];
            is := {};
            if 1 < nops(is1) then
                for i to nops(is1)-1 do
                    ss := is1[i];
                    for j from i+1 to nops(is1) do
                        ss := `charsets/movefactor`(ss,is1[j],ord)
                    od;
                    is := {op(is),`charsets/pfactor`(ss)}
                od;
                {op(is),is1[nops(is1)]}
            else {op(is1)}
            fi
        else {op(is)}
        fi
    end:

# the sequence of distinct factors of f
`charsets/pfactor` :=

    proc(f)
    local i;
        if type(f,integer) then op({})
        elif type(f,`^`) then op(1,f)
        elif type(f,`*`) then
            seq(`charsets/pfactor`(op(i,f)), i = 1 .. nops(f))
        else f
        fi
    end:

# the sequence of factors of f
`charsets/sfactor` :=

    proc(f)
    local i;
        if type(f,integer) then op({})
        elif type(f,`^`) then seq(op(1,f), i = 1 .. op(2,f))
        elif type(f,`*`) then
            seq(`charsets/sfactor`(op(i,f)), i = 1 .. nops(f))
        else f
        fi
    end:

# the set of all nonconstant factors of initials of polys in as: user function
`charsets/iniset` :=

proc(as,ord)
local ind,i;
    ind := 0;
    if nargs <> 2 then ERROR(`wrong number of arguments`)
    elif nops(as) < 1 then ERROR(`no polynomials specified`)
    elif nops(ord) < 1 then ERROR(`no indeterminates specified`)
    elif not type(ord,list) then ERROR(ord,`must be a list`)
    fi;
    if member(false,map(type,ord,name)) then ERROR(`bad variable list`) fi;
    if member(false,map(type,[op(as)],polynom(rational,ord))) then
        ERROR(`input must be polynomials over Q in`,ord)
    fi;
    for i to nops(as) do
        if `charsets/class`(as[i],ord) <= ind then
            ERROR(
    `first argument must be a non-contradictory (weak, quasi-) ascending set`
            )
        else ind := `charsets/class`(as[i],ord)
        fi
    od;
    `charsets/initialset`(expand(as),ord)
end:

# the set of all nonconstant factors of initials of polys in as
`charsets/initialset` :=

    proc(as,ord)
        local i,is,iss;
        is := {};
        for i in as do
            `charsets/initial`(i,ord);
            if `charsets/class`(",ord) > 0 then is := {op(is),"} fi
        od;
        is := `charsets/factorps`(is);
        iss := {};
        for i in is do
            if `charsets/class`(i,ord) > 0 then iss := {op(iss),i} fi
        od;
        iss
    end:

# all irreducible nonconstant factors of a set of polynomials
`charsets/factorps` :=

proc(ps)
local qs,i,j,q;
    qs := {};
    for i in ps do
        q := factor(i);
        if type(q,`*`) then
            for j to nops(q) do
                if not type(op(j,q),integer) then
                    if type(op(j,q),`^`) then
                        qs :=
                            {op(qs),numer(op(1,op(j,q))/lcoeff(op(1,op(j,q))))}
                    else qs := {op(qs),numer(op(j,q)/lcoeff(op(j,q)))}
                    fi
                fi
            od
        elif type(q,`^`) then qs := {op(qs),numer(op(1,q)/lcoeff(op(1,q)))}
        else if not type(q,integer) then qs := {op(qs),numer(q/lcoeff(q))} fi
        fi
    od;
    qs
end:

##### Part II. Routines for Various Decompositions and Others #####

# the char series of polyset ps: user function
`charsets/charser` :=

    proc(ps,lst,medset,y)
    local mset,ord;
        if nargs < 2 then ERROR(`too few arguments`)
        elif nops(ps) < 1 then ERROR(`no polynomials specified`)
        elif nops(lst) < 1 then ERROR(`no indeterminates specified`)
        elif 4 < nargs then ERROR(`too many arguments`)
        fi;
        if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi;
        if member(false,map(type,[op(ps)],polynom(rational,lst))) then
            ERROR(`input must be polynomials over Q in`,lst)
        fi;
        if type(lst,list) then ord := lst
        else ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps)
        fi;
        if 3 < nargs then y := ord fi;
        if nargs < 3 then mset := 'charsetn' else mset := medset fi;
        if member(mset,{'wcharsetn','charsetn','wbasset','trisetc','basset'})
       then
            `charsets/charseries`(expand(ps),ord,`charsets/`.mset)
        else
            ERROR(`medial set must be one of ``basset``,``wbasset``,`.
                ```charsetn``,``wcharsetn`` and ``trisetc```)
        fi
    end:

# the char series of polyset ps
`charsets/charseries` :=

proc(ps,ord,medset)
local qs,cs,iss,n,qsi,qhi,csno,ppi,qqi;
    if type(ps[1],{set,list}) then qhi := {op(ps)} else qhi := {{op(ps)}} fi;
    qsi := {};
    csno := 0;
    ppi := {};
    qqi := {};
    for n from 0 while qhi <> {} do
        qhi := sort([op(qhi)],`charsets/nopsord`);
        qs := qhi[1];
        ppi := `charsets/select`(ppi,nops(qs));
        qqi := {op(ppi[2]),op(qqi)};
        if n = 0 then ppi := {} else ppi := {qs,op(ppi[1])} fi;
        cs := `charsets/charseta`([op(qs)],ord,medset);
        if 1 < printlevel then
            csno := csno+1;
            lprint(
                `characteristic set produced`,csno,nops(qhi),nops(qsi),nops(qs)
                );
            print(cs)
        fi;
        if 0 < `charsets/class`(cs[1],ord) then
            iss := `charsets/initialset`(cs,ord);
            if `charsets/simpa`(iss,cs,ord) <> 0 then qsi := {cs,op(qsi)} fi;
            iss := `charsets/adjoin`(iss,qs,qqi)
        else iss := {} fi;
        if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)}
        else qhi := iss
        fi
    od;
    if qsi <> {} then op(`charsets/contract`(qsi,ord,0)) else [] fi
end:

# the char series of polyset ps -- allowing to remove factors
#       user function
`charsets/mcs` :=

    proc(ps,lst,medset,y)
    local mset,ord;
        if nargs < 2 then ERROR(`too few arguments`)
        elif nops(ps) < 1 then ERROR(`no polynomials specified`)
        elif nops(lst) < 1 then ERROR(`no indeterminates specified`)
        elif 4 < nargs then ERROR(`too many arguments`)
        fi;
        if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi;
        if member(false,map(type,[op(ps)],polynom(rational,lst))) then
            ERROR(`input must be polynomials over Q in`,lst)
        fi;
        if type(lst,list) then ord := lst
        else ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps)
        fi;
        if 3 < nargs then y := ord fi;
        if nargs < 3 then mset := 'charsetn' else mset := medset fi;
        if member(mset,{'wcharsetn','charsetn','wbasset','trisetc','basset'})
        then
            `charsets/fcharser`(expand(ps),ord,`charsets/`.mset)
        else
            ERROR(`medial set must be one of ``basset``,``wbasset``,`.
                ```charsetn``,``wcharsetn`` and ``trisetc```)
        fi
    end:

# the char series of polyset ps -- allowing to remove factors
`charsets/fcharser` :=

proc(ps,ord,medset)
local qs,cs,iss,n,qhi,qsi,factorset,csno,ppi,qqi;
    if type(ps[1],{set,list}) then qhi := {op(ps)} else qhi := {{op(ps)}} fi;
    qsi := {};
    csno := 0;
    ppi := {};
    qqi := {};
    for n from 0 while qhi <> {} do
        qhi := sort([op(qhi)],`charsets/nopsord`);
        qs := qhi[1];
        ppi := `charsets/select`(ppi,nops(qs));
        qqi := {op(ppi[2]),op(qqi)};
        if n = 0 then ppi := {} else ppi := {qs,op(ppi[1])} fi;
        if nops(qs)-3 < nops(ord) then
            cs := `charsets/fcharseta`([op(qs)],ord,medset);
            factorset := op(2,cs[2]);
            cs := cs[1]
        else
            `charsets/charseta`([op(qs)],ord,medset);
            cs := `charsets/removecont`(",ord,'factorset')
        fi;
        if 1 < printlevel then
            csno := csno+1;
            lprint(
                `characteristic set produced`,csno,nops(qhi),nops(qsi),nops(qs)
                );
            print(cs)
        fi;
        if 0 < `charsets/class`(cs[1],ord) then
            iss := `charsets/initialset`(cs,ord);
            if `charsets/simpa`(iss,cs,ord) <> 0 then qsi := {cs,op(qsi)} fi;
            iss := iss union `charsets/factorps`(factorset)
        else iss := `charsets/factorps`(factorset) fi;
        iss := `charsets/adjoin`(iss,qs,qqi);
        if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)}
        else qhi := iss
        fi
    od;
    if qsi <> {} then op(`charsets/contract`(qsi,ord,0)) else [] fi
end:

# the extended char series of polyset ps
#      user function
`charsets/ecs` :=

    proc(ps,lst,medset,y)
    local mset,ord;
        if nargs < 2 then ERROR(`too few arguments`)
        elif nops(ps) < 1 then ERROR(`no polynomials specified`)
        elif nops(lst) < 1 then ERROR(`no indeterminates specified`)
        elif 4 < nargs then ERROR(`too many arguments`)
        fi;
        if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi;
        if type(ps[1],list) then
            if member(false,map(type,ps[1],polynom(rational,lst))) then
                ERROR(`input must be polynomials over Q in`,lst)
            fi
        elif member(false,map(type,[op(ps)],polynom(rational,lst))) then
            ERROR(`input must be polynomials over Q in`,lst)
        fi;
        if type(lst,list) then ord := lst
        else
            if type(ps[1],{set,list}) then
                ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps[1])
            else ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps)
            fi
        fi;
        if 3 < nargs then y := ord fi;
        if nargs < 3 then mset := 'charsetn' else mset := medset fi;
        if member(mset,{'wcharsetn','charsetn','wbasset','trisetc','basset'})
        then
            `charsets/excharser`(expand(ps),ord,`charsets/`.mset)
        else
            ERROR(`medial set must be one of ``basset``,``wbasset``,`.
                ```charsetn``,``wcharsetn`` and ``trisetc```)
        fi
    end:

# the extended char series of polyset ps
`charsets/excharser` :=

proc(ps,ord,medset)
local qs,cs,is,iss,i,j,qsi,qhi,r,rr;
    if type(ps[1],{set,list}) then qhi := {ps} else qhi := {[ps,1]} fi;
    qsi := {};
    while qhi <> {} do
        qs := qhi[1][1];
        cs := `charsets/charseta`([op(qs)],ord,medset);
        if 1 < printlevel then
            lprint(`characteristic set produced`); print(cs)
        fi;
        if 0 < `charsets/class`(cs[1],ord) then
            is := `charsets/initialset`(cs,ord);
            rr := `charsets/nopower`(`charsets/prod`({op(is),qhi[1][2]}));
            `charsets/premas`(rr,cs,ord);
            r := `charsets/simp`(",cs,ord);
            if r <> 0 then
                if r = 1 then qsi := {cs,op(qsi)}
                else qsi := {op(qsi),[cs,`charsets/simpb`(r,rr)]}
                fi
            fi
        else is := []
        fi;
        iss := {};
        if nops(ord) <= nops(ps)+1 then
            for i in is do  iss := {op(iss),[{op(qs),i},qhi[1][2]]} od
        else
            for i to nops(is) do
                if i = 1 then 1 else product('is[j]','j' = 1 .. i-1) fi;
                iss := {op(iss),[{op(qs),is[i]},"*qhi[1][2]]}
            od
        fi;
        if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)}
        else qhi := iss
        fi
    od;
    if qsi <> {} then op(qsi) else [] fi
end:

# simplify r to rr so that Zero(cs/r) = Zero(cs/rr) holds still
`charsets/simp` :=

    proc(r,cs,ord)
    local rr,i,fs,j,ind;
        if r = 0 then 0
        else
            fs := {`charsets/pfactor`(r)};
            rr := 1;
            for i in fs do
                if `charsets/class`(i,ord) > 0 then
                    ind := 0;
                    for j in cs do
                        subs(i=0,j);
                        if " = 0 then
                            if `charsets/class`(`charsets/movefactor`(j,i,ord),
                                ord) = 0 then ind := -1; break fi
                        elif `charsets/class`(",ord) = 0 then
                            ind := 1; break
                        fi
                    od;
                    if ind = 0 then rr := rr*i
                    elif ind = -1 then break
                    fi
                fi
             od;
             if ind = -1 then 0 else rr fi
          fi
      end:

# check whether Zero(cs/fs) is empty
`charsets/simpa` :=

    proc(fs,cs,ord)
    local i,j,ds;
        if nops(cs) = 1 then 1
        else
            ds := [seq(cs[i], i = 1 .. nops(cs)-1)];
            for i in fs do
                for j in ds do
                    if subs(i=0,j) = 0 then
                        if `charsets/class`(`charsets/movefactor`(j,i,ord),
                            ord) = 0 then RETURN(0)
                        fi
                    fi
                od;
            od;
         1
         fi
     end:

# the simpler one of a and b
`charsets/simpb` :=

    proc(a,b)
        if `charsets/measure`(a) < `charsets/measure`(b) then a
        else b
        fi
    end:

# the measure of complexity of a according to number of terms
`charsets/measure` :=

    proc(a)
    local i;
        if type(a,`^`) then nops(op(1,a))
        elif type(a, `*`) then
            sum('`charsets/measure`(op(i,a))','i' = 1 .. nops(a))
        else nops(a)
        fi
    end:


# the extended char series of polyset ps -- allowing to remove factors
#       user function
`charsets/mecs` :=

    proc(ps,lst,medset,y)
    local mset,ord;
        if nargs < 2 then ERROR(`too few arguments`)
        elif nops(ps) < 1 then ERROR(`no polynomials specified`)
        elif nops(lst) < 1 then ERROR(`no indeterminates specified`)
        elif 4 < nargs then ERROR(`too many arguments`)
        fi;
        if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi;
        if type(ps[1],list) then
            if member(false,map(type,ps[1],polynom(rational,lst))) then
                ERROR(`input must be polynomials over Q in`,lst)
            fi
        elif member(false,map(type,[op(ps)],polynom(rational,lst))) then
            ERROR(`input must be polynomials over Q in`,lst)
        fi;
        if type(lst,list) then ord := lst
        else
            if type(ps[1],{set,list}) then
                ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps[1])
            else ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps)
            fi
        fi;
        if 3 < nargs then y := ord fi;
        if nargs < 3 then mset := 'charsetn' else mset := medset fi;
        if member(mset,{'wcharsetn','charsetn','wbasset','trisetc','basset'})
        then
            `charsets/fexcharser`(expand(ps),ord,`charsets/`.mset)
        else
            ERROR(`medial set must be one of ``basset``,``wbasset``,`.
                ```charsetn``,``wcharsetn`` and ``trisetc```)
        fi
    end:

# the extended char series of polyset ps -- allowing to remove factors
`charsets/fexcharser` :=

proc(ps,ord,medset)
local qs,cs,is,iss,n,i,j,qhi,qsi,r,rr,factorset;
    if type(ps[1],{set,list}) then qhi := {ps} else qhi := {[ps,1]} fi;
    qsi := {};
    for n from 0 while qhi <> {} do
        qs := qhi[1][1];
        if n = 0 then
            cs := `charsets/fcharseta`([op(qs)],ord,medset);
            factorset := op(2,cs[2]);
            cs := cs[1]
        else
            `charsets/charseta`([op(qs)],ord,medset);
            cs := `charsets/removecont`(",ord,'factorset')
        fi;
        if 1 < printlevel then
            lprint(`characteristic set produced`); print(cs)
        fi;
        if 0 < `charsets/class`(cs[1],ord) then
            is :=
             `charsets/initialset`(cs,ord) union `charsets/factorps`(factorset)
             ;
            rr := `charsets/nopower`(`charsets/prod`({op(is),qhi[1][2]}));
            `charsets/premas`(rr,cs,ord);
            r := `charsets/simp`(",cs,ord);
            if r <> 0 then
                if r = 1 then qsi := {cs,op(qsi)}
                else qsi := {[cs,`charsets/simpb`(r,rr)],op(qsi)}
                fi
            fi
        else is := `charsets/factorps`(factorset)
        fi;
        iss := {};
        if nops(ord) <= nops(ps)+1 then
            for i in is do  iss := {op(iss),[{op(qs),i},qhi[1][2]]} od
        else
            for i to nops(is) do
                if i = 1 then 1 else product('is[j]','j' = 1 .. i-1) fi;
                iss := {op(iss),[{op(qs),is[i]},"*qhi[1][2]]}
            od
        fi;
        if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)}
        else qhi := iss
        fi
    od;
    if qsi <> {} then op(qsi) else [] fi
end:

# the irreducible char series of polyset ps: user function
`charsets/ics` :=

proc(ps,lst,medset,y)
local mset,ord;
    if nargs < 2 then ERROR(`too few arguments`)
    elif nops(ps) < 1 then ERROR(`no polynomials specified`)
    elif nops(lst) < 1 then ERROR(`no indeterminates specified`)
    elif 4 < nargs then ERROR(`too many arguments`)
    fi;
    if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi;
    if member(false,map(type,[op(ps)],polynom(rational,lst))) then
        ERROR(`input must be polynomials over Q in`,lst)
    fi;
    if type(lst,list) then ord := lst
    else ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps)
    fi;
    if 3 < nargs then y := ord fi;
    if nargs < 3 then mset := 'charsetn' else mset := medset fi;
    if member(mset,{'charsetn','trisetc','basset'}) then
        `charsets/irrcharser`(expand(ps),ord,`charsets/`.mset)
    else
ERROR(`medial set must be one of ``basset``,``charsetn```.` and ``trisetc```)
    fi
end:

# the irreducible char series of polyset ps
#     using new factorization method if m=1, Hu-Wang's method if m=-1
#     and normalized char set if m=0
`charsets/irrcharser` :=

proc(ps,ord,medset,m)
local qs,cs,is,iss,n,i,ts,qsi,qhi,pi,factorset,ppi,qqi,csno,ind,fset,var,mind;
options remember;
    if nargs = 3 then ind := 1 else ind := m fi;
    if type(ps[1],{set,list}) then qhi := {op(ps)} else qhi := {ps} fi;
    qsi := {};
    pi := {};
    csno := 0;
    ppi := {};
    qqi := {};
    var := {seq(op(indets(ps[i])), i = 1 .. nops(ps))};
    if medset = `charsets/basset` then mind := true else mind := false fi;
    for n from 0 while qhi <> {} do
        qhi := sort([op(qhi)],`charsets/nopsord`);
        qs := qhi[1];
        ppi := `charsets/select`(ppi,nops(qs));
        qqi := {op(qqi),op(ppi[2])};
        if n = 0 then ppi := {} else ppi := {op(ppi[1]),qs} fi;
        if nops(qs)-3 < nops(ord) then
            if not mind then
                cs := `charsets/f`.(substring(medset,10 .. length(medset)))(
                    qs,ord,[{},var],'fset');
                if nops(indets(cs[1])) > 2 then
                    cs := `charsets/removecont`(cs,ord,'factorset');
                    factorset := factorset union fset[1]
                else factorset := fset[1]
                fi
            else
                cs := `charsets/fcharseta`(qs,ord,medset);
                factorset := op(2,cs[2]);
                cs := cs[1];
                if 1 < printlevel and ind = 1 then
                    csno := csno+1;
                    lprint(`characteristic set produced`,csno,nops(qhi),
                        nops(qsi),nops(qs));
                    print(cs)
                fi
            fi;
            if ind = 0 then
                cs := [`charsets/fcnormal`(cs,ord)];
                if 1 < nops(cs) then factorset := factorset union op(2,cs[2])
                fi;
                cs := cs[1]
            fi;
            if mind and nops(indets(cs[1])) > 2 then
                cs := `charsets/removecont`(cs,ord,'ts');
                factorset := factorset union ts
            fi
        elif mind then
            cs := `charsets/removecont`(
                `charsets/charseta`([op(qs)],ord,medset),ord,'factorset')
        else
            cs := `charsets/removecont`(medset([op(qs)],ord),ord,'factorset')
        fi;
        if 0 < `charsets/class`(cs[1],ord) then
            ts := `charsets/irras`(cs,ord,ind);
            if ts[2] = 0 then
                if not mind then
                    cs := `charsets/charseta`({op(cs),op(qs)},ord,medset);
                    if 1 < printlevel and ind = 1 then
                        csno := csno+1;
                        lprint(`characteristic set produced`,csno,nops(qhi),
                            nops(qsi),nops(qs));
                        print(cs)
                    fi
                fi;
                if not member(cs,pi) then
                    pi := {cs,op(pi)};
                    if 0 < `charsets/class`(cs[1],ord) then
                        ts := `charsets/irras`(cs,ord,ind);
                        if ts[2] = 0 then
                            qsi := {cs,op(qsi)};
                            if nops(cs) = nops(ord) then
                                is := `charsets/factorps`(factorset)
                            else
                                is := `charsets/initialset`(cs,ord) union
                                    `charsets/factorps`(factorset)
                            fi;
                            iss := `charsets/adjoin`(is,qs,qqi)
                        fi
                    else
                        iss := `charsets/adjoin`(
                            `charsets/factorps`(factorset),qs,qqi)
                    fi
                else
                    iss :=
                       `charsets/adjoin`(`charsets/factorps`(factorset),qs,qqi)
                fi
            fi;
            if ts[2] <> 0 then
                is := `charsets/factorps`(factorset) union ts[1];
                if 1 < ts[2] then
                    is :=
                      is union `charsets/initialset`([op(1 .. ts[2]-1,cs)],ord)
                fi;
                iss := `charsets/adjoin`(is,qs,qqi)
            fi
        else iss := `charsets/adjoin`(`charsets/factorps`(factorset),qs,qqi)
        fi;
        if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)}
        else qhi := iss
        fi
    od;
    if qsi <> {} then op(`charsets/contract`(qsi,ord,1)) else [] fi
end:

# examine the irreducibility of as for irrcharser
`charsets/irras` :=

  proc(as,ord,inda)
  local ind,i,j,qs,n,fs,ja;
  options remember,system;
      ind := 1;
      ja := 0;
      for i to nops(as) do
          fs := `charsets/factorps`({as[i]});
          qs := {};
          for j to nops(fs) do
              if 0 < `charsets/class`(fs[j],ord) then qs := {op(qs),fs[j]} fi
          od;
          if degree(qs[1],`charsets/lvar`(as[i],ord)) <
              degree(as[i],`charsets/lvar`(as[i],ord)) then
              ja := 1; ind := 0; break
          fi
      od;
      if ind = 1 and 1 < nops(as) then
          for n while
              n < nops(as) and degree(as[n],`charsets/lvar`(as[n],ord)) = 1 do

          od;
          if n < nops(as) then qs := `charsets/irrassub`(as,n,ord,inda,'ja')
          else ja := 0
          fi
      fi;
      [qs,ja]
  end:

# subroutine for irras
`charsets/irrassub` :=

 proc(as,n,ord,ind,ja)
 local m,qs,i;
     for m from n+1 while
         m <= nops(as) and degree(as[m],`charsets/lvar`(as[m],ord)) = 1 do

     od;
     if m <= nops(as) then
         if ind = -1 then
             qs := `charsets/factoras`([seq(as[i], i = 1 .. m-1)],as[m],ord)
         else
             qs :=
                 `charsets/newfactoras`([seq(as[i], i = 1 .. m-1)],as[m],ord)
         fi;
         if qs = as[m] then qs := `charsets/irrassub`(as,m,ord,ind,'ja')
         else qs := `charsets/factorps`({qs}); ja := m
         fi
     else ja := 0
     fi;
     qs
 end:

# normalize ascending set cs wrt ord
`charsets/fcnormal` :=

proc(cs,ord)
local n,ini,i,j,ggg,gg,ff,ccs,dd,cd,fs,nt;
    n := nops(cs);
    if n < 2 then cs
    else
        dd := cs[n];
        nt := nops(expand(dd));
        for i from n-1 by -1 to 1 do
            ini := `charsets/initial`(dd,ord);
            if 0 < degree(ini,`charsets/lvar`(cs[i],ord)) then
                ggg := `charsets/gcdex`(cs[i],ini,`charsets/lvar`(cs[i],ord));
                gg := ggg[3];
                if 0 < degree(gg,`charsets/lvar`(cs[i],ord)) then
                    ff := cs[i];
                    gg := {`charsets/pfactor`(gg)};
                    cd := {};
                    for j to nops(gg) do
                        if `charsets/class`(gg[j],ord) =
                            `charsets/class`(cs[i],ord) then
                            ff := `charsets/nopower`(
                                `charsets/movefactor`(ff,gg[j],ord));
                            cd := {op(cd),gg[j]}
                        fi
                    od;
                    if `charsets/class`(ff,ord) = 0 then ccs := [[1]]
                    else
                        ccs := [`charsets/fcnormal`(subs(cs[i] = ff,cs),ord)]
                    fi;
                    if nops(ccs) = 1 then
                        RETURN(ccs[1],`common divisors` = cd)
                    else
                        RETURN(
                           ccs[1],`common divisors` = {op(cd),op(op(2,ccs[2]))}
                           )
                    fi
                else
                    dd :=
                    `charsets/prem`(dd*ggg[2],cs[i],`charsets/lvar`(cs[i],ord))
                    ;
                    dd :=
                    `charsets/movefactor`(dd,`charsets/initial`(cs[i],ord),ord)
                    ;
                    dd := `charsets/nopower`(dd);
                    if 3*nt < nops(expand(dd)) then RETURN(cs) fi
                fi
            fi
        od;
        ccs := [seq(cs[i], i = 1 .. n-1)];
        fs := {`charsets/pfactor`(content(dd,`charsets/lvar`(dd,ord),'dd'))};
        gg := {};
        for i to nops(fs) do
            if 0 < `charsets/class`(fs[i],ord) then gg := {fs[i],op(gg)} fi
        od;
        gg := `charsets/prod`(gg);
        ini := `charsets/initialset`(ccs,ord);
        for i to nops(ini) do  gg := `charsets/movefactor`(gg,ini[i],ord) od;
        dd := `charsets/nopower`(gg)*dd;
        gg := [`charsets/pfactor`(numer(dd/lcoeff(dd)))];
        dd := 1;
        for i to nops(gg) do
            if 0 < `charsets/class`(gg[i],ord) then dd := dd*gg[i] fi
        od;
        if 2*nt < nops(expand(dd)) then cs else [op(ccs),dd] fi
    fi
end:

# the modified gcdex for fcnormal
`charsets/gcdex` := proc(A,B,x)
                  local m,pm,cc,cd,c,c1,c2,d,d1,d2,r,r1,r2,q,II,g;
                  options remember,system;
                      if A = 0 then RETURN([0,1,B]) fi;
                      if B = 0 then RETURN([1,0,A]) fi;
                      cc := content(A,x,c);
                      cd := content(B,x,d);
                      II := readlib(`gcd/degrees`)(c,d,{x},'c','d');
                      pm := 1;
                      c1 := 1;
                      c2 := 0;
                      d1 := 0;
                      d2 := 1;
                      while d <> 0 do
                          r := prem(c,d,x,'m','q');
                          divide(r,pm,'r');
                          divide(m*c1-q*d1,pm,'r1');
                          divide(m*c2-q*d2,pm,'r2');
                          c := d;
                          c1 := d1;
                          c2 := d2;
                          d := r;
                          d1 := r1;
                          d2 := r2;
                          pm := m
                      od;
                      lcoeff(g);
                      subs(II,[c1*cd/",c2*cc/",c*cc*cd/"])
                  end:


# the extended irreducible char series of polyset ps: user function
`charsets/eics` :=

proc(ps,lst,medset,y)
local mset,ord;
    if nargs < 2 then ERROR(`too few arguments`)
    elif nops(ps) < 1 then ERROR(`no polynomials specified`)
    elif nops(lst) < 1 then ERROR(`no indeterminates specified`)
    elif 4 < nargs then ERROR(`too many arguments`)
    fi;
    if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi;
    if type(ps[1],list) then
        if member(false,map(type,ps[1],polynom(rational,lst))) then
            ERROR(`input must be polynomials over Q in`,lst)
        fi
    elif member(false,map(type,[op(ps)],polynom(rational,lst))) then
        ERROR(`input must be polynomials over Q in`,lst)
    fi;
    if type(lst,list) then ord := lst
    else
        if type(ps[1],{set,list}) then
            ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps[1])
        else ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps)
        fi
    fi;
    if 3 < nargs then y := ord fi;
    if nargs < 3 then mset := 'charsetn' else mset := medset fi;
    if member(mset,{'charsetn','trisetc','basset'}) then
        `charsets/exirrcharser`(expand(ps),ord,`charsets/`.mset)
    else
ERROR(`medial set must be one of ``basset``,`.```charsetn`` and ``trisetc```)
    fi
end:

# the extended irreducible char series of polyset ps
`charsets/exirrcharser` :=

proc(ps,ord,medset)
local qs,cs,is,iss,n,i,j,qhi,qsi,r,rr,factorset,mind,var,fset,ind,ts;
    if type(ps[1],{set,list}) then qhi := {ps} else qhi := {[ps,1]} fi;
    if medset = `charsets/basset` then mind := true else mind := false fi;
    qsi := {};
    var := {seq(op(indets(ps[i])), i = 1 .. nops(ps))};
    for n from 0 while qhi <> {} do
        qs := qhi[1][1];
        if not mind then
            if n < 20 then
                `charsets/f`.(substring(medset,10 .. length(medset)))(
                    qs,ord,[{},var],'fset');
                cs:=`charsets/removecont`(",ord,'factorset');
                factorset := factorset union fset[1]
            else
                `charsets/`.(substring(medset,10 .. length(medset)))(
                    qs,ord);
                cs:=`charsets/removecont`(",ord,'factorset')
            fi
        else
            if n < 20 then
                cs := `charsets/fcharseta`([op(qs)],ord,medset);
                factorset := op(2,cs[2]);
                cs := `charsets/removecont`(cs[1],ord,'ts');
                factorset := factorset union ts
            else
                `charsets/charseta`([op(qs)],ord,medset);
                cs:=`charsets/removecont`(",ord,'factorset')
            fi;
            if 1 < printlevel then
                lprint(`characteristic set produced`); print(cs)
            fi
        fi;
        if 0 < `charsets/class`(cs[1],ord) then
            ts := `charsets/irras`(cs,ord,ind);
            if ts[2] = 0 then
                if not mind then
                    cs := `charsets/charseta`({op(cs),op(qs)},ord,medset);
                    if 1 < printlevel then
                        lprint(`characteristic set produced`); print(cs)
                    fi
                fi;
                if 0 < `charsets/class`(cs[1],ord) then
                    ts := `charsets/irras`(cs,ord,ind);
                    if ts[2] = 0 then
                       is := `charsets/initialset`(cs,ord) union
                       `charsets/factorps`(factorset);
                       if nops(cs)=nops(ord) then
                           rr := `charsets/nopower`(qhi[1][2])
                       else rr := `charsets/nopower`(
                           `charsets/prod`({op(is),qhi[1][2]})) fi;
                       `charsets/premas`(rr,cs,ord);
                       r := `charsets/simp`(",cs,ord);
                       if r <> 0 then
                           if r =1 then qsi := {cs,op(qsi)}
                           else qsi := {[cs,`charsets/simpb`(r,rr)],op(qsi)}
                           fi
                       fi
                    fi
                else is := `charsets/factorps`(factorset); ts := [1,0]
                fi
            fi;
            if ts[2] <> 0 then
                if 1 < ts[2] then
                    is := `charsets/initialset`({op(1 .. ts[2]-1,cs)},ord)
                        union `charsets/factorps`(factorset)
                else is := `charsets/factorps`(factorset) union ts[1]
                fi
            fi
        else is := `charsets/factorps`(factorset); ts := [1,0]
        fi;
        iss := {};
        if nops(ord) <= nops(ps)+1 then
            for i in is do  iss := {op(iss),[{op(qs),i},qhi[1][2]]} od
        else
            for i to nops(is) do
                if i = 1 then 1 else product('is[j]','j' = 1 .. i-1) fi;
                iss := {op(iss),[{op(qs),is[i]},"*qhi[1][2]]}
            od
        fi;
        if 1 < ts[2] and ts[1] <> {} then
            for i in ts[1] do
                iss :=
                    {[{op(cs),i},`charsets/prod`({op(is),qhi[1][2]})],op(iss)}
            od
        fi;
        if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)}
        else qhi := iss
        fi
    od;
    if qsi <> {} then op(qsi) else [] fi
end:

# subroutine for irrcharser, qirrcharser and others
`charsets/select` :=

    proc(ppi,n)
    local i,pp,qq;
        pp := {};
        qq := {};
        for i in ppi do
            if n <= nops(i) then qq := {op(qq),i} else pp := {op(pp),i} fi
        od;
        [pp,qq]
    end:

# subroutine for irrcharser, qirrcharser and others
`charsets/adjoin` := proc(is,qs,qh)
                   local iss,i,j,ind,qhi,itt;
                       iss := {};
                       qhi := qh minus {qs};
                       if is <> {} then
                           for i in is do
                               itt := {op(qs),i};
                               ind := 0;
                               if 0 < nops(qhi) then
                                   for j in qhi while ind = 0 do
                                       if itt union j = itt then ind := 1 fi
                                   od
                               fi;
                               if ind = 0 then iss := {op(iss),itt} fi
                           od
                       fi;
                       iss
                   end:

# subroutine for trisersub
`charsets/adjoina` := proc(is,qs,qh)
                    local iss,i,j,ind,qhi,itt;
                        iss := {};
                        qhi := qh minus {qs};
                        if is <> {} then
                            for i in is do
                                itt := {op(qs),i};
                                ind := 0;
                                if 0 < nops(qhi) then
                                    for j in qhi while ind = 0 do
                                        if itt union j = itt then ind := 1 fi
                                    od
                                fi;
                                if ind = 0 then iss := {op(iss),[i,op(qs)]} fi
                            od
                        fi;
                        iss
                    end:

# subroutine for irrcharser, qirrcharser and others
`charsets/nopsord` := proc(a,b)
                    options remember,system;
                        if nops(b) < nops(a) then true else false fi
                    end:

# remove some redundant ascending sets in cs
#     irr=1 for irrcharser, irr=2 for qirrcharser, irr=-1 for trisersub
#     and irr=0 for others
`charsets/contract` :=

   proc(cs,ord,irr)
   local i,j,mem,ts;
       mem := {};
       ts := {};
       if nops(cs) < 2 then cs
       else
           for i to nops(cs)-1 do
               if not member(i,mem) then
                   for j from i+1 to nops(cs) do
                       if not member(j,mem) then
                           if `charsets/linas`(cs[i],ord,irr) and
                               `charsets/contractsub`(cs[i],cs[j],ord) then
                               ts := {cs[j],op(ts)}; mem := {op(mem),j}
                           else
                               if `charsets/linas`(cs[j],ord,irr) and
                                   `charsets/contractsub`(cs[j],cs[i],ord) then
                                   ts := {cs[i],op(ts)}
                               fi
                           fi
                       fi
                   od
               fi
           od;
           {op(cs)} minus ts
       fi
   end:

# check whether all polys in cs1 have remainders 0 wrt cs2
#      but none of their initials does: subroutine for contract
`charsets/contractsub` :=

    proc(cs1,cs2,ord)
    local i,is;
        for i in cs1 do
            if `charsets/premas`(i,cs2,ord) <> 0 then RETURN(false) fi
        od;
        is := `charsets/initialset1`(cs1,ord);
        for i in is do
            if `charsets/premas`(i,cs2,ord) = 0 then RETURN(false) fi
        od;
        true
    end:

# check the irreducibility of as: subroutine for contract
`charsets/linas` :=

proc(as,ord,irr)
local i,j,n,m;
    if irr = 1 then true
    elif irr = 2 then
        if 1 < nops(as) then
            for n while
                n < nops(as) and degree(as[n],`charsets/lvar`(as[n],ord)) = 1
                do

            od;
            if n < nops(as) then
                for m from n+1 while
                 m <= nops(as) and degree(as[m],`charsets/lvar`(as[m],ord)) = 1
                  do

                od;
                if m <= nops(as) then RETURN(false) fi
            fi
        fi;
        true
    else
        for i in as do
            if 1 < degree(i,`charsets/lvar`(i,ord)) then RETURN(false) fi
        od;
        if irr=0 or nops(as)<2 then true
        else
           for i from 2 to nops(as) do
               for j from 1 to i do
               if has(`charsets/initial`(as[i],ord),`charsets/lvar`(as[j],ord))
                   then RETURN(false) fi
               od
           od;
           true
        fi
    fi
end:

# the quasi-irreducible char series of polyset ps: user function
`charsets/qics` :=

    proc(ps,lst,medset,y)
    local mset,ord;
        if nargs < 2 then ERROR(`too few arguments`)
        elif nops(ps) < 1 then ERROR(`no polynomials specified`)
        elif nops(lst) < 1 then ERROR(`no indeterminates specified`)
        elif 4 < nargs then ERROR(`too many arguments`)
        fi;
        if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi;
        if member(false,map(type,[op(ps)],polynom(rational,lst))) then
            ERROR(`input must be polynomials over Q in`,lst)
        fi;
        if type(lst,list) then ord := lst
        else ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps)
        fi;
        if 3 < nargs then y := ord fi;
        if nargs < 3 then mset := 'charsetn' else mset := medset fi;
        if member(mset,{'wcharsetn','charsetn','wbasset','trisetc','basset'})
        then
            `charsets/qirrcharser`(expand(ps),ord,`charsets/`.mset)
        else
            ERROR(`medial set must be one of ``basset``,``wbasset``,`.
                ```charsetn``,``wcharsetn`` and ``trisetc```)
        fi
    end:

# the quasi-irreducible char series of polyset ps
`charsets/qirrcharser` :=

proc(ps,ord,medset)
local qs,cs,is,iss,n,i,ts,qsi,qhi,pi,factorset,ppi,qqi,csno,var,fset,mind;
    qhi := {ps};
    qsi := {};
    pi := {};
    csno := 0;
    ppi := {};
    qqi := {};
    var := {seq(op(indets(ps[i])), i = 1 .. nops(ps))};
    if medset <> `charsets/basset` and medset <> `charsets/wbasset` then
        mind := true
    else mind := false
    fi;
    for n from 0 while qhi <> {} do
        qhi := sort([op(qhi)],`charsets/nopsord`);
        qs := qhi[1];
        ppi := `charsets/select`(ppi,nops(qs));
        qqi := {op(qqi),op(ppi[2])};
        if n = 0 then ppi := {} else ppi := {op(ppi[1]),qs} fi;
        if nops(qs)-3 < nops(ord) then
            if mind then
                cs := `charsets/f`.(substring(medset,10 .. length(medset)))(
                    qs,ord,[{},var],'fset');
                if nops(indets(cs[1])) > 2 then
                    cs := `charsets/removecont`(cs,ord,'factorset');
                    factorset := factorset union fset[1]
                else factorset := fset[1]
                fi
            else
                cs := `charsets/fcharseta`(qs,ord,medset);
                factorset := op(2,cs[2]);
                if 1 < printlevel then
                    csno := csno+1;
                    lprint(`characteristic set produced`,csno,nops(qhi),
                        nops(qsi),nops(qs));
                    print(cs[1])
                fi;
                cs := `charsets/removecont`(cs[1],ord,'ts');
                factorset := factorset union ts
            fi
        elif (not mind) and nops(indets(cs[1])) > 2 then
            cs := `charsets/removecont`(
                `charsets/charseta`([op(qs)],ord,medset),ord,'factorset')
        else
            cs := `charsets/removecont`(medset([op(qs)],ord),ord,'factorset')
        fi;
        if 0 < `charsets/class`(cs[1],ord) then
            ts := `charsets/qirras`(cs,ord);
            if ts[2] = 0 then
                if mind then
                    cs := `charsets/charseta`({op(cs),op(qs)},ord,medset);
                    if 1 < printlevel then
                        csno := csno+1;
                        lprint(`characteristic set produced`,csno,nops(qhi),
                            nops(qsi),nops(qs));
                        print(cs)
                    fi
                fi;
                if not member(cs,pi) then
                    pi := {cs,op(pi)};
                    if 0 < `charsets/class`(cs[1],ord) then
                        ts := `charsets/qirras`(cs,ord);
                        if ts[2] = 0 then
                            qsi := {cs,op(qsi)};
                            if nops(cs) = nops(ord) then
                                is := `charsets/factorps`(factorset)
                            else
                                is := `charsets/initialset`(cs,ord) union
                                    `charsets/factorps`(factorset)
                            fi;
                            iss := `charsets/adjoin`(is,qs,qqi)
                        fi
                    else
                        iss := `charsets/adjoin`(
                            `charsets/factorps`(factorset),qs,qqi)
                    fi
                else
                    iss :=
                       `charsets/adjoin`(`charsets/factorps`(factorset),qs,qqi)
                fi
            fi;
            if ts[2] <> 0 then
                is := `charsets/factorps`(factorset) union ts[1];
                if 1 < ts[2] then
                    is :=
                      is union `charsets/initialset`([op(1 .. ts[2]-1,cs)],ord)
                fi;
                iss := `charsets/adjoin`(is,qs,qqi)
            fi
        else iss := `charsets/adjoin`(`charsets/factorps`(factorset),qs,qqi)
        fi;
        if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)}
        else qhi := iss
        fi
    od;
    if qsi <> {} then op(`charsets/contract`(qsi,ord,2)) else [] fi
end:

# examine the irreducibility of as for qirrcharser
`charsets/qirras` :=

proc(as,ord)
local ind,i,qs,n,m,ja;
options remember,system;
    qs := [];
    ind := 1;
    ja := 0;
    for i to nops(as) do
        qs := `charsets/factorps`({as[i]});
        if numer(qs[1]/lcoeff(qs[1])) <> numer(as[i]/lcoeff(as[i])) then
            ja := 1; ind := 0; break
        fi
    od;
    if ind = 1 and 1 < nops(as) then
        for n while
            n < nops(as) and degree(as[n],`charsets/lvar`(as[n],ord)) = 1 do

        od;
        if n < nops(as) then
            for m from n+1 while
                m <= nops(as) and degree(as[m],`charsets/lvar`(as[m],ord)) = 1
                do

            od;
            if m <= nops(as) then
                lprint(
                 `Warning: factorization over algebraic field required for ics`
                 )
            fi
        fi
    fi;
    [qs,ja]
end:

# factorize poly f over algebraic number filed with minimal polys in as
#       wrt ord: user function
`charsets/cfactor` :=

proc(f,as,ord)
local ind,i,ff;
    if nargs = 1 then RETURN(factor(f)) fi;
    ind := 0;
    ff := numer(f);
    if nargs = 2 then ERROR(`inproper number of arguments`)
    elif nops(as) < 1 then ERROR(`no polynomials specified`)
    elif nops(ord) < 1 then ERROR(`no indeterminates specified`)
    elif not type(ord,list) then ERROR(ord,`must be a list`)
    elif 4 < nargs then ERROR(`too many arguments`)
    fi;
    if member(false,map(type,ord,name)) then ERROR(`bad variable list`) fi;
    if member(false,map(type,[op(as)],polynom(rational,ord))) then
        ERROR(`input must be polynomials over Q in`,ord)
    fi;
    for i to nops(as) do
        if `charsets/class`(as[i],ord) <= ind then
            ERROR(`second argument must be a non-contradictory ascending set`)
        else ind := `charsets/class`(as[i],ord)
        fi
    od;
    lprint(`Warning: Be sure the ascending set is irreducible`);
    if `charsets/class`(ff,ord) <= `charsets/class`(as[nops(as)],ord) then
        factor(f)
    else
        ind := [`charsets/sfactor`(factor(ff))];
        `charsets/@con` := [1961713*_wdm+7131961*_dmw+1371961*_mwd+19901025,
            7131961*_dmw+1371961*_mwd+19901025];
        ind := `charsets/prod`(
            [seq(`charsets/newfactoras`(as,ind[i],ord), i = 1 .. nops(ind))])
            ;
        lcoeff(ff,`charsets/lvar`(ff,ord))*ind/
            lcoeff(ind,`charsets/lvar`(ind,ord));
        if type(",`*`) then "/factor(denom(f)) else f fi
    fi
end:

# factorize poly f over algebraic number filed with minimal polys in as
#       -- new method in combination with Hu-Wang's
`charsets/@con` := [1961713*_wdm+7131961*_dmw+1371961*_mwd+19901025,
                    7131961*_dmw+1371961*_mwd+19901025]:
`charsets/newfactoras` :=

proc(as,f,ord)
local aas,con,i;
    if `charsets/class`(f,ord) = 0 then RETURN(f)
    elif degree(f,`charsets/lvar`(f,ord)) = 1 then RETURN(f)
    fi;
    aas := [];
    con := 2;
    for i to nops(as) do
        degree(as[i],`charsets/lvar`(as[i],ord));
        if 1 < " then con := con,"; aas := [op(aas),expand(as[i])] fi
    od;
    if {con,degree(f,`charsets/lvar`(f,ord))} = {2} and nops(aas) < 3 then
        `charsets/factoras`(aas,f,ord)
    else
        if 1 < printlevel then
            lprint(`newfactoras: factorization over algebraic field -- degree `
                ,degree(f,`charsets/lvar`(f,ord)))
        fi;
        `charsets/newfactorassub`(aas,f,ord)
    fi
end:


# main subroutine for newfactoras
`charsets/newfactorassub` :=

proc(as,f,ord)
local nord,mord,cs,cr,i,j,con,ff,ind,nas,xf,der,ct,fset,fs;
options remember;
    nas := nops(as);
    xf := `charsets/lvar`(f,ord);
    if 1 < nas then
        for i from 2 to nas do
            `charsets/newfactorassub`([seq(as[j], j = 1 .. i-1)],as[i],ord)
                 := as[i]
        od
    fi;
    nord := [xf,seq(`charsets/lvar`(as[i],ord), i = 1 .. nas)];
    mord := [seq(`charsets/lvar`(as[i],ord), i = 1 .. nas),xf];
    con := 0;
    for i from 2 to nops(nord) do  con := con+degree(f,nord[i]) od;
    if con = 0 then
        for i to nas do
            if 1 < degree(as[i],nord[i+1]) then con := con,as[i] fi
        od;
        if nops([con]) = 2 and degree(con[2],`charsets/lvar`(con[2],ord)) = 2
             and degree(f,nord[1]) = 3 then
            RETURN(f)
        fi;
        der := sum('(i-1)*nord[i]','i' = 2 .. nops(nord));
        ff := expand(subs(nord[1] = der+`charsets/@zz`,f));
        nord := subs(nord[1] = `charsets/@zz`,nord);
        con := 1
    elif divide(f,`charsets/@con`[2]) then
        der := sum('`charsets/die`()*nord[i]','i' = 2 .. nops(nord));
        ff := expand(subs(nord[1] = der+`charsets/@zz`,f));
        nord := subs(nord[1] = `charsets/@zz`,nord);
        con := 1
    else ff := f
    fi;
    cs := `charsets/fqcharsetn`([op(as),ff],nord,[{},{}],'fs');
    cs := [`charsets/sfactor`(factor(cs[1])),op(fs[1])];
    if nops(cs) = 1 then
        if 0 < `charsets/class`(cs[1],nord) then
            cs := [`charsets/irrcharser`(
                [op(as),ff,op(cs)],nord,`charsets/charsetn`,-1)]
        else
            cs := [`charsets/irrcharser`(
                [op(as),ff,cs[1]],nord,`charsets/charsetn`,-1)]
        fi;
        ind := 1
    elif divide(f,`charsets/@con`[2]) and divide(f,`charsets/@con`[1]) then
        cs := [`charsets/irrcharser`(
            {seq([op(as),ff,cs[i]], i = 1 .. nops(cs))},nord,
            `charsets/charsetn`,-1)];
        ind := 1
    else ind := 0
    fi;
    `charsets/@con` := [`charsets/@con`[2],f];
    if con = 1 then
        cs := expand(subs(`charsets/@zz` = xf-der,cs))
    fi;
    cr := [];
    ct := [];
    for i to nops(cs) do
        if 1 < printlevel then
            lprint(`characteristic set:`,nops(cs),i)
        fi;

        if ind = 0 then
            [`charsets/fcnormal`(`charsets/fcharsetn`([op(as),f,cs[i]],mord,
            [{},{op(mord)}],'fset'),mord)];
            if nops("[1]) = nas+1 then
                "[1][nas+1];
                if degree(",xf)=1 then ct := [op(ct),"]
                else cr := [op(cr),"]
                fi
            fi
        else
            `charsets/fcharsetn`([op(as),f,op(cs[i])],mord,[{},{op(mord)}],
            'fset');
            if nops(") = nas+1 then ct := [op(ct),"[nas+1]] fi
        fi
    od;
    fset := false;
    [op(cr),op(ct)];
    if " <> [] then
        `charsets/prod`([seq("[i], i = 1 .. nops("))]);
        if `charsets/premas`(
             lcoeff(",`charsets/lvar`(",ord))*f-lcoeff(f,xf)*",as,ord) = 0
             then fset := true
        fi;
    fi;
    if fset then
        if nops(ct) = 0 then
            `charsets/prod`([seq(`charsets/newfactorassub`(as,cr[i],ord),
            i = 1 .. nops(cr))])
        elif nops(cr) = 0 then
            `charsets/prod`([seq(ct[i], i = 1 .. nops(ct))])
        else
            `charsets/prod`([seq(ct[i], i = 1 .. nops(ct))])*
            `charsets/prod`([seq(`charsets/newfactorassub`(as,cr[i],ord),
            i = 1 .. nops(cr))])
        fi
    else
        cr := [];
        for i to nops(cs) do
            if ind = 0 then
                if 1 < printlevel then
                    lprint(`irreducible characteristic series:`,nops(cs),i)
                fi;
                `charsets/irrcharser`([op(as),f,cs[i]],mord,
                `charsets/charsetn`,0)
            else
                if 1 < printlevel then
                    lprint(`characteristic set:`,nops(cs),i)
                fi;
                op(1,[`charsets/fcharseta`(
                    [op(as),f,op(cs[i])],mord,`charsets/`.charsetn)])
            fi;
            if " <> [] then cr := [op(cr),"] fi
        od;
        cr := `charsets/prod`([seq(cr[i][nas+1], i = 1 .. nops(cr))]);
        if `charsets/premas`(
            lcoeff(cr,`charsets/lvar`(cr,ord))*f-lcoeff(f,xf)*cr,as,ord) = 0
            then cr
        else `charsets/factoras`(args)
        fi
    fi
end:

# generate random numbers between 1 and 4
`charsets/die` := rand(1..4):

# factorize poly f over algebraic number filed with minimal polys in as
#       -- Hu-Wang's method
`charsets/factoras` :=

  proc(pas,pf,ord)
  local df,r,s,ff,i,t,fact,gg,hh,fg,z,ind,as,f,sol;
      ind := 0;
      as := [];
      r := 0;
      f := expand(pf);
      for i to nops(pas) do
          df := degree(pas[i],`charsets/lvar`(pas[i],ord));
          if 1 < df then r := r+1; `charsets/@m`.r := df; as := [op(as),pas[i]]
          fi
      od;
      z := `charsets/lvar`(f,ord);
      df := degree(f,z);
      if df = 1 then f
      elif r = 0 then f
      else
          if 1 < printlevel then
              lprint(`factoras: factorization over algebraic field -- degree `,
                  degree(f,`charsets/lvar`(f,ord)))
          fi;
          for s to trunc(1/2*df) while ind = 0 do
              for i to s do
                  `charsets/g`.i := `charsets/summ`(
                      `charsets/@g`[i,seq(`charsets/@k`.t, t = 1 .. r)]*
                      product('`charsets/lvar`(as[t],ord)^`charsets/@k`.t',
                              't' = 1 .. r),r)
              od;
              for i to df-s do
                  `charsets/h`.i := `charsets/summ`(
                      `charsets/@h`[i,seq(`charsets/@k`.t, t = 1 .. r)]*
                      product('`charsets/lvar`(as[t],ord)^`charsets/@k`.t',
                              't' = 1 .. r),r)
              od;
              `charsets/g`.0 := 1;
              `charsets/h`.0 := 1;
              gg := sum('`charsets/g`.i*z^(s-i)','i' = 0 .. s);
              hh := sum('`charsets/h`.i*z^(df-s-i)','i' = 0 .. df-s);
              ff := f-lcoeff(expand(f),z)*expand(gg*hh);
              ff := `charsets/premas`(ff,as,ord);
              fact := {};
              for i from 0 to df-1 do
                  fact :=
                      {op(fact),op(`charsets/coeff`(as,{coeff(ff,z,i)},r,ord))}
              od;
              sol := [`charsets/solveps`(fact,`charsets/getvars`(fact))];
              fg := f;
              if sol <> [] then
                  fg := subs(sol[1],gg)*
                      `charsets/factoras`(as,numer(subs(sol[1],hh)),ord);
                  ind := 1
              fi
          od;
          numer(fg)
      fi
  end:

# subroutine for factoras
`charsets/summ` :=

    proc(ss,r)
        if r = 1 then sum(ss,`charsets/@k`.r = 0 .. `charsets/@m`.r-1)
        else sum(`charsets/summ`(ss,r-1),`charsets/@k`.r = 0 ..
                 `charsets/@m`.r-1)
        fi
    end:

# subroutine for factoras
`charsets/coeff` :=

    proc(as,ss,r,ord)
    local k,i,j,qs;
        qs := ss;
        for j from r by -1 to 1 do
            qs := {seq(seq(coeff(qs[i],`charsets/lvar`(as[j],ord),k),
                k = 0 .. `charsets/@m`.j-1), i = 1 .. nops(qs))}
        od;
        qs minus {0}
    end:

# subroutine for factoras
`charsets/getvars` :=

    proc(as)
    local ind,ind1,i;
        if type(as,{set,list}) then
            {seq(op(`charsets/getvars`(as[i])), i = 1 .. nops(as))}
        else
            ind := {};
            ind1 := indets(as);
            for i in ind1 do
                if type(i,indexed) then
                    if op(0,i) = `charsets/@g` or op(0,i) = `charsets/@h` then
                        ind := {op(ind),i}
                    fi
                fi
            od;
            ind
        fi
    end:


# find rational zeros of polyset ps
`charsets/solveps` :=

   proc(ps,lst)
   local cs,ord,sol,j,phi,qs,qs1,n,factorset;
   options remember;
       if 1 < printlevel then
           lprint(`solveps: trying rational solutions of equations`,
               op(`charsets/index`([op(ps)],[op(lst)])))
       fi;
       ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps);
       sol := {};
       cs := `charsets/fqcharsetn`(ps,ord,[{},{}],'factorset');
       factorset := factorset[1];
       phi := {ps};
       for n while phi <> {} do
           if sol <> {} then break fi;
           if 1 < n then
               cs := `charsets/charseta`(phi[1],ord,`charsets/`.wcharsetn);
               factorset := {}
           fi;
           sol := {op(sol),`charsets/solveasr`(cs,ord,'qs1')};
           if n = 1 then sol := `charsets/verify`(sol,ps,ord) fi;
           qs := `charsets/factorps`(qs1) union `charsets/factorps`(factorset);
           if qs <> {} then
               if 1 < nops(phi) then
                   phi := {op(2 .. nops(phi),phi),
                       seq([op(phi[1]),qs[j]], j = 1 .. nops(qs))}
               else phi := {seq([op(phi[1]),qs[j]], j = 1 .. nops(qs))}
               fi
           else
               if 1 < nops(phi) then phi := {op(2 .. nops(phi),phi)}
               else phi := {}
               fi
           fi
       od;
       if sol = {} then op({}) else op(sol) fi
   end:

# subroutine for solveps
`charsets/verify` :=

    proc(sol,ps,ord)
    local i,j,sss;
        for i to nops(sol) do
            if simplify(subs(sol[i],ps)) = {0} then RETURN({sol[i]})
            else
                {seq(op(1,sol[i][j])-op(2,sol[i][j]), j = 1 .. nops(sol[i]))};
                sss := {`charsets/solveps`(ps union ",ord)};
                if sss <> {} then RETURN(sss) fi
            fi
        od;
        {}
    end:

# prepare a list of triangular forms from polyset ps: user function
`charsets/triser` :=

proc(ps,lst,y)
local i,ord;
    if nargs < 1 then ERROR(`too few arguments`)
    elif nops(ps) < 1 then ERROR(`no polynomials specified`)
    elif 3 < nargs then ERROR(`too many arguments`)
    elif nargs = 2 then
        if nops(lst) < 1 then ERROR(`no indeterminates specified`) fi
    fi;
    if nargs = 2 then
        if member(false,map(type,lst,name)) then ERROR(`bad variable list`)
        fi
    fi;
    if member(false,map(type,[op(ps)],polynom(rational,lst))) then
        ERROR(`input must be polynomials over Q in`,lst)
    fi;
    if nargs < 2 then
        ord := `charsets/reorder`(
            [seq(op(indets(ps[i])), i = 1 .. nops(ps))],`charsets/degord`,ps)
    elif type(lst,list) then ord := lst
    else ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps)
    fi;
    if 2 < nargs then y := ord fi;
    `charsets/trisersub`(expand(ps),ord)
end:

# prepare a list of triangular forms from polyset ps
`charsets/trisersub` :=

proc(ps,ord)
local qs,cs,iss,n,i,qhi,qsi,factorset,csno,ppi,qqi,ind,mem;
options remember;
    ind := 0;
    for i to nops(ps) do
        if nops(expand(ps[i])) < 3 then ind := 1; break fi
    od;
    if ind = 1 then
        cs := `charsets/fcharseta`([op(ps)],ord,`charsets/`.charsetn)
    else cs := `charsets/fcharseta`([op(ps)],ord,`charsets/`.qcharsetn)
    fi;
    factorset := op(2,cs[2]);
    cs := cs[1];
    qhi := {{op(ps)}};
    qsi := {};
    csno := 0;
    ppi := {};
    qqi := {};
    for n from 0 while qhi <> {} do
        qhi := sort([op(qhi)],`charsets/nopsord`);
        qs := qhi[1];
        ppi := `charsets/select`(ppi,nops(qs));
        qqi := {op(qqi),op(ppi[2])};
        if n = 0 then ppi := {}
        else
            ppi := {op(ppi[1]),qs};
            ind := 0;
            for i to nops(ps) do
                if nops(expand(ps[i])) < 3 then ind := 1; break fi
            od;
            if ind = 1 then
                cs := `charsets/nopower`(
                    `charsets/charseta`(qs,ord,`charsets/`.charsetn));
                factorset := {}
            elif qs <> mem and 4 < degree(qs[1],ord) then
                cs := `charsets/fcharseta`(qs,ord,`charsets/`.qcharsetn);
                factorset := op(2,cs[2]);
                cs := cs[1]
            elif nops(qs)-3 < nops(ord) then
                cs := `charsets/fcharseta`(qs,ord,`charsets/`.wcharsetn);
                factorset := op(2,cs[2]);
                cs := cs[1]
            else
                cs := `charsets/nopower`(
                    `charsets/charseta`(qs,ord,`charsets/`.wcharsetn));
                factorset := {}
            fi
        fi;
        mem := qs;
        if 1 < printlevel then
            csno := csno+1;
            lprint(
                `characteristic set produced`,csno,nops(qhi),nops(qsi),nops(qs)
                );
            print(cs)
        fi;
        if 0 < `charsets/class`(cs[1],ord) then
            iss := `charsets/initialset`(cs,ord);
            if `charsets/simpa`(iss,cs,ord) <> 0 then qsi := {cs,op(qsi)} fi;
            iss := iss union `charsets/factorps`(factorset)
        else iss := `charsets/factorps`(factorset) fi;
        iss := `charsets/adjoina`(iss,qs,qqi);
        if 1 < nops(qhi) then qhi := {op(iss),op(2 .. nops(qhi),qhi)}
        else qhi := iss
        fi
    od;
    if qsi <> {} then op(`charsets/contract`(qsi,ord,-1)) else {} fi
end:

# solve a set of poly eqs ps=0: user function
`charsets/csolve` :=

proc(ps,lst,y)
local i,ord,qsi;
    if nargs < 1 then ERROR(`too few arguments`)
    elif nops(ps) < 1 then ERROR(`no polynomials specified`)
    elif 3 < nargs then ERROR(`too many arguments`)
    elif nargs = 2 then
        if nops(lst) < 1 then ERROR(`no indeterminates specified`) fi
    fi;
    if nargs = 2 then
        if member(false,map(type,lst,name)) then ERROR(`bad variable list`)
        fi
    fi;
    if member(false,map(type,[op(ps)],polynom(rational,lst))) then
        ERROR(`input must be polynomials over Q in`,lst)
    fi;
    if nargs < 2 then
        ord := `charsets/reorder`(
            [seq(op(indets(ps[i])), i = 1 .. nops(ps))],`charsets/degord`,ps)
    elif type(lst,list) then ord := lst
    else ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps)
    fi;
    if 2 < nargs then y := ord fi;
    qsi := {`charsets/trisersub`(expand(ps),ord)};
    if qsi = {{}} then {}
    else op({seq(`charsets/solveas`(qsi[i],ord), i = 1 .. nops(qsi))})
    fi
end:

# find zeros of ascending set as
`charsets/solveas` :=

proc(cs,ord)
local is,ss,sol,solm,i,j,k;
    sol := {solve({cs[1]},{`charsets/lvar`(cs[1],ord)})};
    if 1 < nops(cs) then
        for i from 2 to nops(cs) do
            is := `charsets/initial`(cs[i],ord);
            solm := {};
            for j to nops(sol) do
                ss :=
                    {solve({subs(sol[j],cs[i])},{`charsets/lvar`(cs[i],ord)})};
                for k to nops(ss) do
                    if subs(op(sol[j]),ss[k],is) <> 0 then
                        solm := {op(solm),{op(sol[j]),op(ss[k])}}
                    fi
                od
            od;
            sol := solm
        od
    fi;
    op(sol)
end:

# find rational zeros of ascending set cs
`charsets/solveasr` :=

proc(cs,ord,qs)
local is,ss,ts,sol,solm,i,j,k;
    ts := {};
    if 0 < `charsets/class`(cs[1],ord) then
        sol := {`charsets/solvel`(cs[1],`charsets/lvar`(cs[1],ord))};
        if 1 < nops(cs) then
            for i from 2 to nops(cs) do
                if 1 <= nops(sol) then
                    is := `charsets/initial`(cs[i],ord);
                    solm := {};
                    for j to nops(sol) do
                        if subs(op(sol[j]),is) = 0 then ts := {op(ts),is}
                        else
                            ss := {`charsets/solvel`(
                                subs(sol[j],cs[i]),`charsets/lvar`(cs[i],ord))}
                                ;
                            solm := {op(solm),
                               seq({op(sol[j]),op(ss[k])}, k = 1 .. nops(ss))
                               }
                        fi
                    od;
                    sol := solm
                else break
                fi
            od
        fi
    else sol := {}
    fi;
    if 2 < nargs then qs := ts fi;
    op(sol)
end:

# find rational zeros of polynomial f wrt x: subroutine for solveasr
`charsets/solvel` :=

    proc(f,x)
    local g,i,sol;
        sol := {};
        if nops(indets(f)) = 1 then
            g := `charsets/getfactor`(f,x);
            for i in g do  sol := {op(sol),solve({i},{x})} od
        else
            g := `charsets/factorps`({numer(f)});
            for i in g do
                if degree(i,x) = 1 then sol := {op(sol),solve({i},{x})} fi
            od
        fi;
        op(sol)
    end:

# find a list of distinct linear factors of univariate poly f
`charsets/getfactor` :=

  proc(f,x)
  local q,qs,j;
      q := `charsets/getfact`(f,x);
      qs := {};
      if type(q,`*`) then
          for j to nops(q) do
              if not type(op(j,q),integer) then
                  if type(op(j,q),`^`) then
                      qs := {op(qs),numer(op(1,op(j,q))/lcoeff(op(1,op(j,q))))}
                  else qs := {op(qs),numer(op(j,q)/lcoeff(op(j,q)))}
                  fi
              fi
          od
      elif type(q,`^`) then qs := {op(qs),numer(op(1,q)/lcoeff(op(1,q)))}
      else if not type(q,integer) then qs := {op(qs),numer(q/lcoeff(q))} fi
      fi;
      [op(qs)]
  end:

# find the product of linear factors of univar poly f using `factor/linfacts`
`charsets/getfact` :=

    proc(ff,x)
    local i,f;
        if degree(ff,x) = 1 then RETURN(ff) fi;
        f := convert(ff,`sqrfree/sqrfree`,x);
        if type(f,`^`) then
            readlib(factor);
            readlib(`factor/polynom`);
            readlib(`factor/unifactor`);
            readlib(`factor/linfacts`)(expand(op(1,f)),x)
        elif type(f,`*`) then
            {seq(`charsets/getfact`(op(i,f),x), i = 1 .. nops(f))};
            product(op(i,"),i = 1 .. nops("))
        else
            readlib(factor);
            readlib(`factor/polynom`);
            readlib(`factor/unifactor`);
            readlib(`factor/linfacts`)(numer(f),x)
        fi
    end:

# the irreducible decomposition of algebraic variety defined by ps
#      user function
`charsets/ivd` :=

proc(ps,lst,medset,y)
local mset,ord;
    if nargs < 2 then ERROR(`too few arguments`)
    elif nops(ps) < 1 then ERROR(`no polynomials specified`)
    elif nops(lst) < 1 then ERROR(`no indeterminates specified`)
    elif 4 < nargs then ERROR(`too many arguments`)
    fi;
    if member(false,map(type,lst,name)) then ERROR(`bad variable list`) fi;
    if member(false,map(type,[op(ps)],polynom(rational,lst))) then
        ERROR(`input must be polynomials over Q in`,lst)
    fi;
    if type(lst,list) then ord := lst
    else ord := `charsets/reorder`([op(lst)],`charsets/degord`,ps)
    fi;
    if 3 < nargs then y := ord fi;
    if nargs < 3 then mset := 'charsetn' else mset := medset fi;
    if member(mset,{'charsetn','trisetc','basset'}) then
        `charsets/irrvardec`(expand(ps),ord,`charsets/`.mset)
    else
ERROR(`medial set must be one of ``basset``,`.```charsetn`` and ``trisetc```)
    fi
end:

# the irreducible decomposition of algebraic variety defined by ps
`charsets/irrvardec` :=

proc(ps,ord,medset)
local phi,psi,qs,gb,zz,i,j,ts,mem,is,qq;
    qq := nops(ps);
    ts := {};
    mem := {};
    if 1 < printlevel and nargs < 3 then lprint(`variable order chosen:`,ord)
    fi;
    psi := [`charsets/irrcharser`(ps,ord,medset)];
    phi := [];
    for i to nops(psi) do
        if nops(psi[i]) <= qq then phi := [op(phi),psi[i]] fi
    od;
    psi := [];
    if phi <> [[]] then
        for i to nops(phi) do
            if nops(phi[i]) = nops(ord) then is := {}
            else is := `charsets/initialset`(phi[i],ord)
            fi;
            if is <> {} then
                qs :=
                 [op(phi[i]),seq(`charsets/@z`.j*is[j]-1, j = 1 .. nops(is))]
                 ;
                zz := [seq(`charsets/@z`.(nops(is)-j+1), j = 1 .. nops(is))];
                gb := grobner['gbasis'](
                 qs,[op(zz),seq(ord[nops(ord)-j+1], j = 1 .. nops(ord))],
                 'plex');
                qs := [];
                for j to nops(gb) do
                    if {op(zz)} minus indets(gb[j]) = {op(zz)} then
                        qs := [gb[j],op(qs)]
                    fi
                od
            else qs := phi[i]
            fi;
            psi := [op(psi),qs]
        od;
        if 1 < nops(psi) then
            for i to nops(psi)-1 do
                if not member(i,mem) then
                    for j from i+1 to nops(psi) do
                        if not member(j,mem) then
                            if `charsets/remseta`(psi[i],phi[j],ord) = {} then
                                ts := {op(ts),psi[j]}; mem := {j,op(mem)}
                            else
                                if `charsets/remseta`(psi[j],phi[i],ord) = {}
                                     then
                                    ts := {op(ts),psi[i]}
                                fi
                            fi
                        fi
                    od
                fi
            od;
            op({op(psi)} minus ts)
        else psi[1]
        fi
    else []
    fi
end:

# The help files for the individual functions are in the file "help"
`help/text/charsets` := TEXT(
`HELP FOR : The characteristic sets package`,
`   `,
`SYNTAX : charsets[ <function> ]( args );   or   <function> ( args );`,
`   `,
`SYNOPSIS :   `,
`   `,
`- To use a charsets function, type "with( charsets, <function> );" to define`,
`  <function>, or type "with( charsets );"  to define all charsets functions,`,
`  or the longer package format charsets[ <function> ].`,
`   `,
`- The functions available are:`,
`   `,
`        cfactor      charser      charset     csolve      ecs`,
`        eics         mcharset     mcs         mecs        ics`,
`        iniset       ivd          qics        remset      triser`,
`   `,
`- For help with a particular function do:  ?<function>`,
`   `,
`- For example, to compute the characteristic series of a set of polynomials`,
`  PS over the variables  X  =  [x[1], x[2] ,..., x[n]], with respect to the`,
`  ordering x[1] < x[2] < ... < x[n], do either:`,
`   `,
`                  with(charsets, charser);  charser( PS, X );`,
`  or:   `,
`                        charsets[charser]( PS, X ).`,
`   `,
`SEE ALSO : with`
):

# Help pages for the charsets package
`help/text/cfactor` := TEXT(
`FUNCTION: cfactor - factorize polynomial over algebraic number field`,
`   `,
`CALLING SEQUENCE: cfactor(F,AS,X);`,
`   `,
`PARAMETERS:  F          - (multivariate) polynomial`,
`             AS         - irreducible ascending set with respect to X`,
`             X          - list of names`,
`   `,
`SYNOPSIS :   `,
`   `,
`- Let the irreducible ascending set AS and the polynomial F be of the form`,
`  AS  = [A[1](u[1] ,..., u[d],x[1]) ,..., A[r](u[1] ,..., u[d], y[1] ,...,`,
`  y[r])]  and  F = F(u[1] ,..., u[d], y[1] ,..., y[r], y).  cfactor(F,AS,X)`,
`  computes  the factorization of F over the algebraic number  field  Q(u[1]`,
`  ,..., u[d], y[1] ,..., y[r]),  where  Q  denotes the rational field, u[1]`,
`  ,..., u[d] are transcendental elements and y[1] ,..., y[r] are algebraic`,
`  elements with  y[i] being an extended zero of A[i] in  Q(u[1] ,..., u[d],`,
`  y[1] ,..., y[i-1]) for each i.`,
`- cfactor(F,AS,X) - F has its remainder 0 with respect to AS.`,
`- X must be the list of leading variables of polynomials in AS and F, i.e.,`,
`  X := [y[1], y[2]  ,...,  y[r], y], with respect to which [A[1] ,..., A[r],`,
`  F] forms an ascending set.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`AS := [-2*x1*x2+2*x2**2+2*x2**2*x1+1+x1]:`,
`F := -x2*x1**3-3*x1**2*x2-x3**2-x1*x3**2+x3*x1*x2+x2*x1**2*x3:`,
`cfactor(F,AS,[x1,x2,x3]);   yields   (1+x1)*(-x3-x1+2*x1*x2)*(x3-x1+x1*x2)`,
`AS := [-1+b+6*b**2+12*b**3]:`,
`F := -6156+29484*b+23328*c+11664*c**3-17712*c**2+67392*c**2*b**2`,
`     +58320*c**2*b-58320*b*c-50544*b**2*c+17496*b**2:`,
`cfactor(F,AS,[b,c]);   yields`,
`    144*(-27*b-32*c+72*b*c+48*b**2*c+27*b**2+18+27*c**2)*(12*b**2+7*b+3*c-1)`,
`AS := [-5-2*x**2, x+5*y+3*y**3]:`,
`F := 15*z**2*x+1725*z*y+39*y**2*x*z+3450*y*x-1495*x-7475*y-75*z:`,
`cfactor(F,AS,[x,y,z]);   yields  -3*x*(10*x+5*z+13*y**2)*(-z+46*y*x)`,
`   `,
`SEE ALSO : Factor`
):
`help/text/charser` := TEXT(
`FUNCTION: charser - compute (weak) characteristic series`,
`   `,
`CALLING SEQUENCE: charser(PS,X);  or  charser(PS,X,medset);`,
`                  or  charser(PS,X,medset,'Y');`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             X          - set or list of names`,
`                          (not including parameters)`,
`             medset     - medial set: 'basset' (basic set),`,
`                          'wbasset' (weak basic set),`,
`                          'charsetn' (nearly characteristic set - default),`,
`                          'wcharsetn' (weak nearly characteristic set),`,
`                          'trisetc' (strong triangularized set)`,
`             Y          - (a name) list of reordered names`,
`   `,
`SYNOPSIS :   `,
`   `,
`- charser(PS,X,medset)  computes  the  (weak)  characteristic series of the`,
`  polynomial set PS with respect to the variables X.  The  resulting series`,
`  CR  is returned as a sequence of lists of polynomials of the form  CS[1],`,
`  ..., CS[e],  where each  CS[i] is an ascending set if medset is 'basset',`,
`  'charsetn' or 'trisetc';  a weak  ascending  set  if  it  is 'wbasset' or`,
`  'wcharsetn'.`,
`- In any case, the zero relation  Zero(PS) = union(Zero(CS[i]/J[i]), i=1..e)`,
`  holds, where  Zero(PS)  denotes the set of common zeros of polynomials in`,
`  PS,  Zero(CS[i]/J[i]) = Zero(CS[i]) minus Zero(J[i]), J[i] is the product`,
`  of initials of polynomials in CS[i].`,
`- A  list  of  variables  X := [x[1], x[2] ,..., x[n]] induces the ordering`,
`  x[1] < x[2] < ... < x[n].`,
`- If the variables are  given as a set  X := {x[1], x[2] ,..., x[n]},  they`,
`  are automatically reordered  to  be "heuristically optimal";  The list of`,
`  reordered variables is assigned to the fourth argument 'Y' if it appears.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`PS := [x1*x3-3*x2-1, -2*x4*x1-x3*x2*x4+2*x3, x2*x4**2-x4*x1+x3]:`,
`charser(PS,[x1,x2,x3,x4]);      yields`,
`[x1, 3*x2+1, x3*(x3-12), x3*(6+x4)], [x2*(18*x2**2*x1**2+12*x1**2*x2+2*x1**2+`,
`27*x2**4+27*x2**3+9*x2**2+x2+36*x2**2*x1+24*x1*x2+4*x1),-x1*x3+3*x2+1,2*x4*x1`,
`**2+3*x2**2*x4-6*x2+x4*x2-2], [x1, 3*x2+1, x3, x4**2], [x1**2, 3*x2+1, x1*x3,`,
`x4**2*x1]   `,
`charser(PS,{x1,x2,x3,x4},wbasset,'Y');  yields`,
`[x4**2*(6*x4*x1**2+2*x4**3*x1**3-3*x1+3*x4*x1**3-6-x4-2*x1*x4**2+x4**2*x1**`,
`2), -3*x4*x1+3*x3+x4**2*x1*x3-x4**2, -x1*x3+3*x2+1], [x4**3-9, 3+x1*x4**2,`,
`6*x3*x4+x4**2*x3+18+3*x3**2, -x1*x3+3*x2+1], [2*x4**2+3, 12*x4*x1**2-12-2*x4`,
`-3*x1**2 , 2*x4*x1-2*x3+x1*x3-1, -x1*x3+3*x2+1], [x4, x3, -x1*x3+3*x2+1]`,
`Y;  yields   [x4, x1, x3, x2]`,
`   `,
`SEE ALSO : mcs, ecs, mecs, iniset`
):
`help/text/charset` := TEXT(
`FUNCTION: charset - compute (weak, quasi-) characteristic set`,
`   `,
`CALLING SEQUENCE: charset(PS,X);  or  charset(PS,X,medset);`,
`                 or  charset(PS,X,medset,'Y');`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             X          - set or list of names`,
`                         (not including parameters)`,
`             medset     - medial set: 'basset' (basic set),`,
`                          'wbasset' (weak basic set),`,
`                          'qbasset' (quasi-basic set),`,
`                          'charsetn' (nearly characteristic set - default),`,
`                          'wcharsetn' (weak nearly characteristic set),`,
`                          'qcharsetn' (quasi-nearly characteristic set),`,
`                          'triset' (triangularized set),`,
`                          'trisetc' (strong triangularized set)`,
`             Y          - (a name) list of reordered names`,
`   `,
`SYNOPSIS :   `,
`   `,
`- charset(PS,X,medset)  computes  the  (weak, quasi-)  characteristic set of`,
`  the polynomial set PS with respect to the variables X.  The  resulting set`,
`  CS  is returned as a list of polynomials.  CS  is  a characteristic set if`,
`  medset is 'basset',  'charsetn'  or  'trisetc';  a weak characteristic set`,
`  if it is  'wbasset'  or  'wcharsetn';  a quasi-characteristic set if it is`,
`  'qbasset', 'qcharsetn' or 'triset'.`,
`- In any case,  the zero relation  Zero(CS/J) -< Zero(PS) -< Zero(CS) holds,`,
`  where  Zero(PS)  denotes  the  set  of  common zeros of polynomials in PS,`,
`  Zero(CS/J) = Zero(CS) minus Zero(J),  J is  the  product  of  initials  of`,
`  polynomials in  CS  and  -< stands for ``be contained in``.`,
`- A list  of  variables  X := [x[1], x[2] ,..., x[n]]  induces  the ordering`,
`  x[1] < x[2] < ... < x[n].`,
`- If the variables are given as a set X := {x[1], x[2] ,..., x[n]}, they are`,
`  are  automatically  reordered  to  be "heuristically optimal"; The list of`,
`  reordered variables is assigned to the fourth argument 'Y' if it appears.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`PS := [x**2-3*x*z+2, x*y+z**3, 5*y**2-3*z**3-7]:`,
`charset(PS,[x,y,z],basset);      yields`,
`[-4803*x**8-81*x**11-486*x**9-972*x**7-648*x**5+5*x**12+60*x**10+800*x**6+`,
`1200*x**4+960*x**2+320, 27*x**4*y+x**6+6*x**4+12*x**2+8, x**2-3*x*z+2]`,
`charset(PS,{x,y,z},wcharsetn);  yields`,
`[-315*z**8-135*z**11+25*z**12+60*z**9+176*z**6+168*z**3+196, -5*z**6+9*x*z**`,
`4+21*x*z-6*z**3-14, 9*z**7+21*z**4+5*z**6*y+6*y*z**3+14*y]`,
`charset(PS,{x,y,z},triset,'Y');      yields`,
`[-315*z**8-135*z**11+25*z**12+60*z**9+176*z**6+168*z**3+196, -5*z**6+9*x*z**`,
`4+21*x*z-6*z**3-14, x*y+z**3]`,
`Y; yields  [z, x, y]`,
`   `,
`SEE ALSO : mcharset, iniset`
):
`help/text/csolve` := TEXT(
`FUNCTION: csolve - solve system of polynomial equations`,
`   `,
`CALLING SEQUENCE: csolve(PS);  or  csolve(PS,X);`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             X          - set or list of names`,
`                         (not including parameters, default is indets(PS))`,
`   `,
`SYNOPSIS :   `,
`   `,
`- csolve(PS,X)  finds  all  solutions of the system of polynomial equations`,
`  PS = 0 with respect to variables in  X  from the list of triangular forms`,
`  prepared basically by 'triser'.`,
`- If X is given as a list  X := [x[1], x[2] ,..., x[n]], the ordering x[1]`,
`  < x[2] < ... < x[n] is used in the characteristic sets computation.`,
`- If X is  given as a set  X := {x[1], x[2] ,..., x[n]}, the variables are`,
`  are  automatically  reordered  to  be "heuristically optimal".`,
`- If X is omitted, the set indets(F) is used as a default.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`PS := [x*y+2*x**2-y+1, 2*y**2-3*x*y-x-2, 2*x**2-3*y+1+2*y**3-3*y**2*x]:`,
`csolve(PS);  yields`,
`{y = 1, x = 0}, {x = -(-29/28-9/28*I*7**(1/2))/(5/4-9/28*I*7**(1/2)),`,
`y = -3/4+3/28*I*7**(1/2)}, {x = -(-29/28+9/28*I*7**(1/2))/(5/4+9/28*I*7**`,
`(1/2)), y = -3/4-3/28*I*7**(1/2)}`,
`csolve(PS,[x,y]);  yields`,
`{x = 1/4+9/28*I*7**(1/2), y=-(-9/28+9/28*I*7**(1/2))/(-3/4+9/28*I*7**(1/2))},`,
`{y = 1, x = 0}, {x = 1/4-9/28*I*7**(1/2), y=-(-9/28-9/28*I*7**(1/2))/`,
`(-3/4-9/28*I*7**(1/2))}`,
`   `,
`SEE ALSO : triser`
):
`help/text/ecs` := TEXT(
`FUNCTION: ecs - compute extended (weak) characteristic series`,
`   `,
`CALLING SEQUENCE: ecs(PS,X);  or  ecs(PS,X,medset);`,
`                  or  ecs(PS,X,medset,'Y'); or ecs([PS,G],X);`,
`                  or  ecs([PS,G],X,medset);`,
`                  or  ecs([PS,G],X,medset,'Y');`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             G          - non-zero (multivariate) polynomial in X`,
`             X          - set or list names`,
`                        (not including parameters)`,
`             medset     - medial set: 'basset' (basic set),`,
`                          'wbasset' (weak basic set),`,
`                          'charsetn' (nearly characteristic set - default),`,
`                          'wcharsetn' (weak nearly characteristic set),`,
`                          'trisetc' (strong triangularized set)`,
`             Y          - (a name) list of reordered names`,
`   `,
`SYNOPSIS :   `,
`   `,
`- ecs(PS,X,medset)  or  ecs([PS,G],X,medset)  computes  the  extended (weak)`,
`  characteristic series of the polynomial set PS  or [PS,G] with  respect to`,
`  the variables X.  The  resulting  series  CR is returned as a sequence  of`,
`  lists of a polynomial set and a polynomial of the form [CS[1],F[1]] ,...,`,
`  [CS[e],F[e]],  where  F[i]  is  a  non-zero  polynomial,  [CS[i],F[i]]  is`,
`  replaced  by  CS[i] while  Zero(CS[i]/F[i]) is examined to be  Zero(CS[i])`,
`  and each CS[i] is an ascending set if medset is  'basset',  'charsetn'  or`,
`  'trisetc'; a weak ascending set if it is 'wbasset' or 'wcharsetn'.`,
`- In any case, the zero relation  Zero(PS) = union(Zero(CS[i]/F[i]), i=1..e)`,
`  or Zero(PS/G) = union(Zero(CS[i]/F[i]), i = 1..e)  holds,  where  Zero(PS)`,
`  denotes the set of common zeros of polynomials  in  PS  and  Zero(PS/G) =`,
`  Zero(PS) minus Zero(G).`,
`- A list of  variables  X := [x[1], x[2] ,..., x[n]]  induces  the  ordering`,
`  x[1] < x[2] < ... < x[n].`,
`- If the variables are given as a set X := {x[1], x[2] ,..., x[n]}, they are`,
`  are  automatically reordered  to  be "heuristically optimal"; The list of`,
`  reordered variables is assigned to the fourth argument 'Y' if it appears.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`PS := [x1*x3-3*x2-1, -2*x4*x1-x3*x2*x4+2*x3, x2*x4**2-x4*x1+x3]:`,
`ecs([PS,x1],[x1,x2,x3,x4]);      yields`,
`[[x2*(18*x2**2*x1**2+12*x1**2*x2+2*x1**2+27*x2**4+27*x2**3+9*x2**2+x2`,
`+36*x1*x2**2+24*x1*x2+4*x1), x1*x3-3*x2-1, 2*x4*x1**2+3*x2**2*x4-6*x2+x2*x4`,
`-2], (2*x1**2+3*x2**2+x2)*x1]`,
`ecs(PS,{x1,x2,x3,x4},trisetc,'Y');  yields`,
`[x4, x3, 3*x2+1], [[x4**2*(2*x1**3*x4**3+6*x1**2*x4-3*x1+3*x1**3*x4`,
`-2*x4**2*x1-6-x4+x1**2*x4**2), -3*x1*x4+3*x3+x4**2*x3*x1-x4**2, x1*x2*x4**2`,
`-x1**2*x4+3*x2+1], x4*(2*x4**2+3)*(x4**2*x1+3)], [x4**3-9, x4**2*x1+3,`,
`6*x4*x3+18+x4**2*x3+3*x3**2, x4*x3+3+9*x2], [[2*x4**2+3, 12*x1**2*x4-2*x4`,
`-12-3*x1**2, 2*x1*x4-2*x3+ x3*x1-1, x2*(4*x1*x4-8*x4-x1+2)], (x1-2)*(4*x4-1)]`,
`Y;  yields  [x4, x1, x3, x2]`,
`   `,
`SEE ALSO : ecs, mcs, mecs`
):
`help/text/eics` := TEXT(
`FUNCTION: eics - compute extended irreducible characteristic series`,
`   `,
`CALLING SEQUENCE: eics(PS,X);  or  eics(PS,X,medset);`,
`                  or  eics(PS,X,medset,'Y');`,
`                  or  eics([PS,G],X);`,
`                  or  eics([PS,G],X,medset);`,
`                  or  eics([PS,G],X,medset,'Y');`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             G          - non-zero (multivariate) polynomial in X`,
`             X          - set or list of names`,
`                          (not including parameters)`,
`             medset     - medial set: 'basset' (basic set),`,
`                          'charsetn' (nearly characteristic set - default),`,
`                          'trisetc' (strong triangularized set)`,
`             Y          - (a name) list of reordered names`,
`SYNOPSIS :   `,
`   `,
`- eics(PS,X,medset)  or   eics([PS,G],X,medset)   computes   the    extended`,
`  irreducible characteristic series of the polynomial set PS  or [PS,G] with`,
`  respect to the variables X.  The  resulting  series  CR  is  returned as a`,
`  sequence of lists of a  polynomial  set  and  a  polynomial  of  the  form`,
`  [CS[1],F[1]] ,..., [CS[e],F[e]], where  CS[i]  is an irreducible ascending`,
`  set, F[i] is a non-zero polynomial and [CS[i],F[i]] is replaced  by  CS[i]`,
`  while Zero(CS[i]/F[i]) is examined to be Zero(CS[i]).`,
`- The zero relation  Zero(PS) = union(Zero(CS[i]/F[i]), i=1..e) or Zero(PS/G)`,
`  = union(Zero(CS[i]/F[i]), i = 1..e) holds,  where Zero(PS) denotes the set`,
`  of common zeros of polynomials in PS, Zero(PS/G)=Zero(PS) minus Zero(G).`,
`- A list of  variables  X := [x[1], x[2] ,..., x[n]]  induces  the  ordering`,
`  x[1] < x[2] < ... < x[n].`,
`- If the variables are given  as  a  set  X := {x[1], x[2] ,..., x[n]}, they`,
`  are automatically reordered  to  be  "heuristically optimal";  The list of`,
`  reordered variables is assigned to the fourth argument 'Y' if it appears.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`PS := [x1*x3-3*x2-1, -2*x4*x1-x3*x2*x4+2*x3, x2*x4**2-x4*x1+x3]:`,
`eics([PS,3*x2+1],[x1,x2,x3,x4]);      yields`,
`[[3*x2**2+x2+2*x1**2+4*x1, x3*x1-3*x2-1, 3*x2+1+2*x1*x4], x1*(x1+2)],`,
`[x2, x3*x1-1, x1**2*x4-1]`,
`eics([PS,x1*x3+3*x3],{x1,x2,x3,x4},basset,'Y');  yields`,
`[[2*x4**2*x1+x4+6+3*x1, x3+2*x4, 2*x2*x4**2-4*x4+3*x2+1], x4*(2*x4**2+3)`,
`*(-9+x4**3)*(x1**2*x4-1)*(x3*x1+3*x3)], [[x1**2*x4-1, x1*x4-x3, x2], 3*x1*x4`,
`+1], [[-9+x4**3, 3*x1+x4, x3+2*x4, -9*x2-3+2*x4**2], 6*x4**2+x4-15]`,
`Y;  yields  [x4, x1, x3, x2]`,
`   `,
`SEE ALSO : ics, qics`
):
`help/text/ics` := TEXT(
`FUNCTION: ics - compute irreducible characteristic series`,
`   `,
`CALLING SEQUENCE: ics(PS,X);  or  ics(PS,X,medset);`,
`                  or  ics(PS,X,medset,'Y');`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             X          - set or list of names`,
`                          (not including parameters)`,
`             medset     - medial set: 'basset' (basic set),`,
`                          'charsetn' (nearly characteristic set - default),`,
`                          'trisetc' (strong triangularized set)`,
`             Y          - (a name) list of reordered names`,
`   `,
`SYNOPSIS :   `,
`   `,
`- ics(PS,X,medset)  computes  the  irreducible characteristic series of the`,
`  polynomial set PS with respect to the variables X.  The  resulting series`,
`  CR is returned as  a  sequence  of  lists  of  polynomials  of  the  form`,
`  CS[1] ,..., CS[e],  where  each  CS[i]  is  an irreducible ascending set.`,
`- The zero relation  Zero(PS) = union(Zero(CS[i]/J[i]), i = 1 .. e)  holds,`,
`  where Zero(PS) denotes the set of common  zeros  of  polynomials  in  PS,`,
`  Zero(CS[i]/J[i]) =  Zero(CS[i]) minus Zero(J[i]),  J[i]  is  the  product`,
`  of initials of polynomials in CS[i].`,
`- A list of  variables  X := [x[1], x[2] ,..., x[n]]  induces  the ordering`,
`  x[1] < x[2] < ... < x[n].`,
`- If the variables are  given  as  a set X := {x[1], x[2] ,..., x[n]}, they`,
`  are  automatically reordered to  be "heuristically optimal"; The list  of`,
`  reordered variables is assigned to the fourth argument 'Y' if it appears.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`PS := [x1*x3-3*x2-1, -2*x4*x1-x3*x2*x4+2*x3, x2*x4**2-x4*x1+x3]:`,
`ics(PS,[x1,x2,x3,x4]);      yields`,
`[3*x2+1, x3, x4], [2*x1**2+4*x1+3*x2**2+x2, -x1*x3+3*x2+1, 3*x2+1+2*x4*x1],`,
`[x1, 3*x2+1, x3-12, 6+x4], [x2, x1*x3-1, x4*x1**2-1]`,
`ics(PS,{x1,x2,x3,x4},basset,'Y');  yields`,
`[x4, x3, 3*x2+1], [2*x1*x4**2+x4+6+3*x1, x3+2*x4, 2*x2*x4**2-4*x4+3*x2+1],`,
`[x4*x1**2-1, x4*x1-x3, x2]`,
`Y;  yields  [x4, x1, x3, x2]`,
`   `,
`SEE ALSO : qics, eics, iniset`
):
`help/text/iniset` := TEXT(
`FUNCTION: iniset - compute initial set`,
`   `,
`CALLING SEQUENCE: iniset(AS,X);`,
`   `,
`PARAMETERS:  AS         - (weak, quasi-) ascending set with respect to X`,
`             X          - list of names`,
`   `,
`SYNOPSIS :   `,
`   `,
`- iniset(AS,X)  computes  the  set  of all distinct factors of initials of`,
`  polynomials in a  (weak, quasi-) ascending set  AS  with  respect to the`,
`  variable ordering X.`,
`- A list of variables  X := [x[1], x[2]  ,...,  x[n]] induces the ordering`,
`  x[1] < x[2] < ... < x[n] .`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`iniset([u*x^3-5*x-2,x^3*y-2*x^2*y-13,y*x*z^3-y*z^2-x],[x,y,z]);   yields`,
`{y, x, x - 2}`,
`AS := [-315*z**8-135*z**11+25*z**12+60*z**9+176*z**6+168*z**3+196, 5*z**6`,
`-9*x*z**4-21*x*z+6*z**3+14, 9*z**7+21*z**4+5*z**6*y+6*y*z**3+14*y]:`,
`iniset(AS,[z,x,y]);      yields     [z, 5*z**6+6*z**3+14, 3*z**3+7]`
):
`help/text/ivd` := TEXT(
`FUNCTION: ivd - compute irreducible variety decomposition`,
`   `,
`CALLING SEQUENCE: ivd(PS,X);  or  ivd(PS,X,medset);`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             X          - set or list of names`,
`                          (not including parameters)`,
`             medset     - medial set: 'basset' (basic set),`,
`                          'charsetn' (nearly characteristic set - default),`,
`                          'trisetc' (strong triangularized set)`,
`   `,
`SYNOPSIS :   `,
`   `,
`- ivd(PS,X,medset)  computes the irreducible decomposition of the algebraic`,
`  variety defined by polynomial set PS with respect to the variables X. The`,
`  result VS is returned as a sequence of polynomial sets of the form VS[1],`,
`  ..., VS[e], where each VS[i] defines an irreducible algebraic variety.`,
`- The zero relation  Zero(PS) = union(Zero(VS[i]), i = 1..e)  holds,  where`,
`  Zero(PS) denotes the set of common zeros of polynomials in  PS.`,
`- If  X is given as a list X := [x[1], x[2] ,..., x[n]],  the ordering x[1]`,
`  < x[2] < ... < x[n] is used in the characteristic set computation.`,
`- If the variables  are  given  as  a  set  X := {x[1], x[2] ,..., x[n]}, a`,
`  "heuristically optimal"  ordering  of  them is used in all characteristic`,
`  sets computation.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`PS := [x1*x3-3*x2-1, -2*x4*x1-x3*x2*x4+2*x3, x2*x4**2-x4*x1+x3]:`,
`ivd(PS,[x1,x2,x3,x4]);      yields`,
`[x2, x1*x3-1, x4-x3**2], [3*x2+1, x3, x4], [3*x2**2+x2+2*x1**2+4*x1,`,
`x1*x3-3*x2-1, 4+2*x1+x3*x2, x3+2*x4]`,
`ivd(PS,{x1,x2,x3,x4},basset);  yields`,
`[x4, x3, 3*x2+1], [2*x1*x4**2+x4+6+3*x1, x3+2*x4, 3*x2+1+2*x4*x1], [x4*x1**2`,
`-1, -x4*x1+x3, x2]`
):
`help/text/mcharset` := TEXT(
`FUNCTION: mcharset - compute modified (weak, quasi-) characteristic set`,
`   `,
`CALLING SEQUENCE: mcharset(PS,X);  or  mcharset(PS,X,medset,'Y');`,
`                  or  mcharset(PS,X,medset,'Y');`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             X          - set or list of names`,
`                        (not including parameters)`,
`             medset     - medial set: 'basset' (basic set),`,
`                          'wbasset' (weak basic set),`,
`                          'qbasset' (quasi-basic set),`,
`                          'charsetn' (nearly characteristic set - default),`,
`                          'wcharsetn' (weak nearly characteristic set),`,
`                          'qcharsetn' (quasi-nearly characteristic set),`,
`                          'triset' (triangularized set),`,
`                          'trisetc' (strong triangularized set)`,
`             Y          - (a name) list of reordered names`,
`   `,
`SYNOPSIS :   `,
`   `,
`- mcharset(PS,X,medset) computes the modified  (weak, quasi-) characteristic`,
`  set  CS  of  the  polynomial  set  PS with respect to the variables X. The`,
`  result is returned as a list of  polynomials  and a set of factors removed`,
`  during the computation of the  form  CS,  factors removed = FS.  CS  is  a`,
`  characteristic set if medset is  'basset', 'charsetn' or 'trisetc'; a weak`,
`  characteristic  set   if  it  is  'wbasset'   or  'wcharsetn';  a  quasi-`,
`  characteristic set if it is  'qbasset',  'qcharsetn' or 'triset'.`,
`- In any case,  the  zero  relations  Zero(CS/J) -< Zero(PS), Zero(PS/F) -<`,
`  Zero(CS)  hold,  where  Zero(PS)  denotes  the  set  of  common  zeros  of`,
`  polynomials in PS, Zero(CS/J) = Zero(CS) minus Zero(J),  J  is the product`,
`  of initials of polynomials in CS,  F is the product of all elements in  FS`,
`  and -< stands for ``be contained in``.`,
`- A list of  variables  X := [x[1], x[2] ,..., x[n]]  induces  the  ordering`,
`  x[1] < x[2] < ... < x[n].`,
`- If the variables are given as  a  set  X := {x[1], x[2] ,..., x[n]}, they`,
`  are  automatically reordered  to  be "heuristically optimal"; The list of`,
`  reordered variables is assigned to the fourth argument 'Y' if it appears.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charset):`,
`PS := [x**2-3*x*z+2, x*y+z**3, 5*y**2-3*z**3]:`,
`mcharset(PS,[z,y,x]);      yields`,
`[-135*z**5+25*z**6+60*z**3+36, 15*y*z+5*z**3+6, -15*z**4+5*x*z**3+6*x],`,
`factors removed = {z}`,
`mcharset(PS,{y,z,x},qbasset,'Y');  yields`,
`[-135*z**5+25*z**6+60*z**3+36, -5*z**3+9*x*z-6, x*y+z**3]`,
`factors removed = {z}`,
`Y;  yields  [z, x, y]`,
`mcharset(PS,[y,x,z],trisetc);      yields`,
`[5832+15625*y**6+151875*y**5+33750*y**4+24300*y**2,5*y+3*x,18+45*y*z+25*y**2]`,
`factors removed = {x, y}`,
`   `,
`SEE ALSO : charset, iniset`
):
`help/text/mcs` := TEXT(
`FUNCTION: mcs - compute (weak) characteristic series`,
`   `,
`CALLING SEQUENCE: mcs(PS,X);  or  mcs(PS,X,medset);`,
`                  or  mcs(PS,X,medset,'Y');`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             X          - set or list of names`,
`                          (not including parameters)`,
`             medset     - medial set: 'basset' (basic set),`,
`                          'wbasset' (weak basic set),`,
`                          'charsetn' (nearly characteristic set - default),`,
`                          'wcharsetn' (weak nearly characteristic set),`,
`                          'trisetc' (strong triangularized set)`,
`             Y          - (a name) list of reordered names`,
`   `,
`SYNOPSIS :   `,
`   `,
`- mcs(PS,X,medset)  computes  the  (weak)  characteristic   series   of  the`,
`  polynomial set PS with respect to the variables X.  The  resulting  series`,
`  CR  is  returned  as  a  sequence  of  lists  of  polynomials  of the form`,
`  CS[1] ,..., CS[e], where each CS[i] is  an  ascending  set  if  medset  is`,
`  'basset', 'charsetn' or 'trisetc'; a weak ascending set if it is 'wbasset'`,
`  or 'wcharsetn'.`,
`- In any case, the zero relation  Zero(PS) = union(Zero(CS[i]/J[i]), i=1..e)`,
`  holds, where  Zero(PS)  denotes the set of common zeros  of polynomials in`,
`  PS, Zero(CS[i]/J[i]) = Zero(CS[i]) minus Zero(J[i]),  J[i]  is the product`,
`  of initials of polynomials in CS[i].`,
`- mcs  has  a  same  functionality  as  that  of  charser.  But  during  the`,
`  computation  of characteristic sets in  mcs,  some  possible  factors  are`,
`  examined and allowed to be removed.`,
`- A list of  variables  X := [x[1], x[2] ,..., x[n]]  induces  the  ordering`,
`  x[1] < x[2] < ... < x[n].`,
`- If the variables are given as  a  set  X := {x[1], x[2] ,..., x[n]},  they`,
`  are  automatically reordered  to  be "heuristically optimal"; The  list of`,
`  reordered variables is assigned to the fourth argument 'Y' if it appears.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`PS := [x1*x3-3*x2-1, -2*x4*x1-x3*x2*x4+2*x3, x2*x4**2-x4*x1+x3]:`,
`mcs(PS,[x1,x2,x3,x4]);      yields`,
`[x1, 3*x2+1, x3,x4], [27*x2**4+27*x2**3+9*x2**2+18*x2**2*x1**2+12*x1**2*x2+x2`,
`+2*x1**2+24*x1*x2+4*x1+36*x2**2*x1,-x1*x3+3*x2+1, 2*x4*x1**2+3*x2**2*x4-6*x2+`,
`x4*x2-2], [x1, 3*x2+1, x3-12, 6+x4], [x2, x1*x3-1, x4*x1**2-1]`,
`mcs(PS,{x1,x2,x3,x4},wcharsetn,'Y');  yields`,
`[x4**3-9, 3+x1*x4**2, 6*x3*x4+x4**2*x3+18+3*x3**2, -x1*x3+3*x2+1],`,
`[2*x4**2+3, 12*x4*x1**2-12-2*x4-3*x1**2, 2*x4*x1-2*x3+x1*x3-1, -x1*x3+3*x2`,
`+1], [x4, x3, -x1*x3+3*x2+1], [6*x4*x1**2+2*x4**3*x1**3-3*x1+3*x4*x1**3-6`,
`-x4-2*x1*x4**2+x4**2*x1**2, -3*x4*x1+3*x3+x4**2*x1*x3-x4**2, -x1*x3+3*x2+1]`,
`Y;  yields  [x4, x1, x3, x2]`,
`   `,
`SEE ALSO : charser, ecs, mecs, iniset`
):
`help/text/mecs` := TEXT(
`FUNCTION: mecs - compute extended (weak) characteristic series`,
`   `,
`CALLING SEQUENCE: mecs(PS,X);  or  mecs(PS,X,medset);`,
`                  or  mecs(PS,X,medset,'Y'); mecs([PS,G],X);`,
`                  or  mecs([PS,G],X,medset);`,
`                  or  mecs([PS,G],X,medset,'Y');`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             G          - non-zero (multivariate) polynomial in X`,
`             X          - set or list of names`,
`                        (not including parameters)`,
`             medset     - medial set: 'basset' (basic set),`,
`                          'wbasset' (weak basic set),`,
`                          'charsetn' (nearly characteristic set - default),`,
`                          'wcharsetn' (weak nearly characteristic set),`,
`                          'trisetc' (strong triangularized set)`,
`             Y          - (a name) list of reordered names`,
`   `,
`SYNOPSIS :   `,
`   `,
`- mecs(PS,X,medset)  or  mecs([PS,G],X,medset) computes the  extended  (weak)`,
`  characteristic series of the polynomial set PS  or  [PS,G] with respect to`,
`  the variables X. The resulting series  CR  is  returned  as  a sequence of`,
`  lists of a polynomial set and a polynomial of the form  [CS[1],F[1]] ,...,`,
`  [CS[e],F[e]], where F[i] is a non-zero polynomial, [CS[i],F[i]] is replaced`,
`  by CS[i] while  Zero(CS[i]/F[i])  is  examined to be  Zero(CS[i])  and each`,
`  CS[i] is an ascending set if medset is  'basset', 'charsetn'  or 'trisetc';`,
`  a weak ascending set if it is 'wbasset' or 'wcharsetn'.`,
`- In any case, the zero  relation  Zero(PS) = union(Zero(CS[i]/F[i]), i=1..e)`,
`  or  Zero(PS) = union(Zero(CS[i]/J[i]), i = 1 . .e)  holds,  where  Zero(PS)`,
`  denotes the set of common zeros of  polynomials  in  PS  and  Zero(PS/G) =`,
`  Zero(PS) minus Zero(G).`,
`- mecs has a same functionality as that of ecs. But during  the  computation`,
`  of  characteristic sets  in  mcs,  some possible factors are  examined and`,
`  allowed to be removed.`,
`- A list of  variables  X := [x[1], x[2] ,..., x[n]]  induces  the  ordering`,
`  x[1] < x[2] < ... < x[n].`,
`- If the variables  are  given as a  set  X := {x[1], x[2] ,..., x[n]}, they`,
`  are automatically  reordered  to  be "heuristically optimal"; The  list of`,
`  reordered variables is assigned to the fourth argument 'Y' if it appears.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`PS := [x1*x3-3*x2-1, -2*x4*x1-x3*x2*x4+2*x3, x2*x4**2-x4*x1+x3]:`,
`mecs([PS,x1],[x1,x2,x3,x4]);      yields`,
`[[27*x2**4+27*x2**3+9*x2**2+18*x2**2*x1**2+12*x1**2*x2+x2+2*x1**2+4*x1`,
`+24*x1*x2+36*x1*x2**2, x3*x1-3*x2-1, 2*x1**2*x4+3*x4*x2**2-6*x2+x4*x2-2],`,
`x1*x2*(3*x2**2+x2+2*x1**2)], [x2, x3*x1-1, x1**2*x4-1]`,
`mecs([PS,x2],{x1,x2,x3,x4},trisetc);  yields`,
`[[6*x1**2*x4+2*x1**3*x4**3-3*x1+3*x1**3*x4-6-x4-2*x4**2*x1+x1**2*x4**2,`,
`-3*x1*x4+3*x3+x4**2*x3*x1-x4**2, 3*x2+1+x4**2*x1*x2-x1**2*x4], x4*(-2*x4**2`,
`-3+2*x4**3*x1**2+3*x1**2*x4)], [x4, x3, 3*x2+1], [[-9+x4**3, 3+x4**2*x1,`,
`6*x3*x4+18+x3*x4**2+3*x3**2, x3*x4+3+9*x2], x4*(2*x4**2+3)*x2]`,
`Y; yields  [x4, x1, x3, x2]`,
`   `,
`SEE ALSO : charser, mcs, ecs`
):
`help/text/qics` := TEXT(
`FUNCTION: qics - compute quasi-irreducible characteristic series`,
`   `,
`CALLING SEQUENCE: qics(PS,X);  or  qics(PS,X,medset);`,
`                  or  qics(PS,X,medset,'Y');`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             X          - set or list of names`,
`                          (not including parameters)`,
`             medset     - medial set: 'basset' (basic set),`,
`                          'wbasset' (weak basic set),`,
`                          'charsetn' (nearly characteristic set - default),`,
`                          'charsetn' (weak nearly characteristic set),`,
`                          'trisetc' (strong triangularized set)`,
`             Y          - (a name) list of reordered names`,
`   `,
`SYNOPSIS :   `,
`   `,
`- qics(PS,X,medset)  computes  the  quasi-irreducible  characteristic series`,
`  of the polynomial set PS with respect to the  variables X.  The  resulting`,
`  series CR is returned as a sequence of lists  of  polynomials  of the form`,
`  CS[1] ,..., CS[e], where each CS[i] is  an  ascending  set  if  medset  is`,
`  'basset', 'charsetn' or 'trisetc'; a weak ascending set if it is 'wbasset'`,
`  or 'wcharsetn' and all polynomials in CS[i] are irreducible.`,
`- The zero  relation  Zero(PS) = union(Zero(CS[i]/J[i]), i = 1 .. e)  holds,`,
`  where  Zero(PS)  denotes the set  of  common zeros of polynomials  in  PS,`,
`  Zero(CS[i]/J[i]) = Zero(CS[i]) minus Zero(J[i]),  J[i]  is  the product of`,
`  initials of polynomials in CS[i].`,
`- A list of  variables  X := [x[1], x[2] ,..., x[n]]  induces  the  ordering`,
`  x[1] < x[2] < ... < x[n].`,
`- If the  variables  are given  as a set X := {x[1], x[2] ,..., x[n]}, they`,
`  are  automatically reordered to  be "heuristically optimal"; The list  of`,
`  reordered variables is assigned to the fourth argument 'Y' if it appears.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`PS := [x4**2+x4**2*x1-x4*x2-x1*x2*x4+x1*x2+3*x2, x4*x1+x3-x1*x2,`,
`x4*x3-2*x2**2-x1*x2-1]:`,
`qics(PS,[x1,x2,x3,x4]);      yields`,
`[x1, 2*x2**2+1, x3, x4**2-x2*x4+3*x2], [1+x1, x2, x3-1, -1+x4], [-2*x1*x2`,
`+2*x2**2+2*x1*x2**2+1+x1, -x1**3*x2-3*x1**2*x2-x3**2-x1*x3**2+x1*x2*x3`,
`+x1**2*x2*x3, x4*x1+x3-x1*x2], [1+x1, x2, x3+1, 1+x4]`,
`qics(PS,{x1,x2,x3,x4},wbasset,'Y');  yields`,
`[x4+x2-1, -2*x1*x2+2*x2**2+2*x1*x2**2+1+x1, x4*x1+x3-x1*x2], [x4-2*x2+1,`,
`-2*x1*x2+2*x2**2+2*x1*x2**2+1+x1, x4*x1+x3-x1*x2]`,
`Y;  yields  [x2, x4, x1, x3]`,
`   `,
`SEE ALSO : ics, eics, iniset`
):
`help/text/remset` := TEXT(
`FUNCTION: remset - compute remainder set`,
`   `,
`CALLING SEQUENCE: remset(PS,AS,X);`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             AS         - (weak, quasi-) ascending set with respect to X`,
`             X          - list of names`,
`   `,
`SYNOPSIS :   `,
`   `,
`- remset(PS,AS,X) computes the set of all non-zero remainders of polynomials`,
`  in PS with respect to a (weak, quasi-) ascending set  AS  and the variable`,
`  ordering X.`,
`- A  list  of variables  X := [x[1], x[2]  ,...,  x[n]] induces the ordering`,
`  x[1] < x[2] < ... < x[n].`,
`- If PS contains only one polynomial F, remset(F,AS,X) may be used.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`PS := [x**2*u-x*y+u*z**2, x*y-2*u**2-x**2-2, z**2-u*x**3]:`,
`AS := [x**2-u*x+2, x*y-u*x-2*u**2, x**2*z**2-4*u*z**2-y**2*z**2-11]:`,
`remset(PS,AS,[x,y,z]);      yields`,
`{11+8*u**5*x-8*u**2*x-4*u*x+4*u**3*x+4*u**4*x-8*u**4-8*u**2-8*u**3,`,
`u*(-12*u**4*x-8*u**5-28*u**2*x-16*u**3*x+4*x-15*u*x+8*u**5*x+30+16*u**2`,
`-16*u**4+32*u)}`,
`remset(x**2+2*y**2+3*z**2,AS,[x,y,z]);        yields`,
`120*u**7*x+32*u**8+20*u**5*x-128*x*u**4-101*u**3*x+164*u*x+48*u**2*x`,
`+76*u**6*x-112*u**6-152*u**5-148+80*u**3-72*u**4-32*u+122*u**2`
):
`help/text/triser` := TEXT(
`FUNCTION: triser - computes triangular series`,
`   `,
`CALLING SEQUENCE: triser(PS);  or  triser(PS,X); or triser(PS,X,Y);`,
`   `,
`PARAMETERS:  PS         - set or list of (multivariate) polynomials in X`,
`             X          - set or list of names`,
`                         (not including parameters, default is indets(PS))`,
`             Y          - (a name) list of reordered names`,
`   `,
`SYNOPSIS :   `,
`   `,
`- triser(PS,X,Y)   computes  a  collection  of  triangular  forms  from  the`,
`  polynomial set PS.  The resulting collection  TR is returned as a sequence`,
`  of lists of polynomials  of the form TS[1] ,..., TS[e],  where each  CS[i]`,
`  is an ascending, weak or quasi-ascending set.`,
`- TR satisfies the zero relation Zero(PS) = union(Zero(TS[i]/J[i]), i=1..e),`,
`  where Zero(PS) denotes the set of  common  zeros  of  polynomials  in  PS,`,
`  Zero(TS[i]/J[i]) = Zero(TS[i]) minus Zero(J[i]) and J[i] is the product of`,
`  initials of polynomials in TS[i].`,
`- If X is given as a list  X := [x[1], x[2] ,..., x[n]],  the  ordering x[1]`,
`  < x[2] < ... < x[n] is used in the characteristic sets computation.`,
`- If X is  given as a set X := {x[1], x[2] ,..., x[n]},  the  variables are`,
`  automatically  reordered  to  be  "heuristically optimal";  The  list  of`,
`  reordered variables is assigned to the third argument 'Y' if it appears.`,
`- If X is omitted, the set indets(F) is used as a default.`,
`   `,
`EXAMPLES :   `,
`   `,
`with(charsets):`,
`PS := [x*y+2*x**2-y+1, 2*y**2-3*x*y-x-2, 2*x**2-3*y+1+2*y**3-3*y**2*x]:`,
`triser(PS);      yields`,
`[-7*y**3+3*y+14*y**4-19*y**2+9, 2*y**2-3*x*y-x-2]`,
`triser(PS,[x,y]);      yields`,
`[x, y-1], [14*x**2-7*x+11, x*y+2*x**2-y+1]`,
`   `,
`SEE ALSO : charser, mcs, ecs, qics`
):

save `charsets.m`;
quit
