share/algebra/charsets/charsets.mac

   1 /* CharSets Version 1.0 (December 1990) */
   2 /* CharSets Version 1.1 (January 1992) for Maple V */
   3 /* CharSets Version 1.2 (January 1994) for Maple V.2 */
   4 /* CharSets Version 2.0 (February 1996) for Maple V.3 */
   5 /*#################################################################### */
   6 /*                                                                   # */
   7 /*                  CHARACTERISTIC SETS PACKAGE                      # */
   8 /*                                                                   # */
   9 /*   Author:  Dongming Wang                                          # */
  10 /*            Laboratoire LEIBNIZ                                    # */
  11 /*            Institut IMAG, 46, avenue Felix Viallet                # */
  12 /*            38031 Grenoble Cedex, France                           # */
  13 /*            E-mail: Dongming.Wang@imag.fr                          # */
  14 /*                                                                   # */
  15 /*   Date:    February 1996                                          # */
  16 /*                                                                   # */
  17 /*   Copyright (C) 1990-1996 by Dongming Wang                        # */
  18 /*                                                                   # */
  19 /*   Copyright Notice:  Permission is granted to use, copy or re-    # */
  20 /*            distribute this package, provided that the title is    # */
  21 /*            retained and the file is not altered.                  # */
  22 /*                                                                   # */
  23 /*#################################################################### */
  24 /*===================================================================# */
  25 /*    This package is implemented for computing characteristic sets  # */
  26 /*  of (multivariate) polynomial sets, decomposing  polynomial sets  # */
  27 /*  into ascending sets and irreducible ascending sets, decomposing  # */
  28 /*  algebraic  varieties into  irreducible components,  decomposing  # */
  29 /*  polynomial   ideals   into   primary   components,  factorizing  # */
  30 /*  polynomials over algebraic extension fields and solving systems  # */
  31 /*  of polynomial equations.  It is  based on the method of charac-  # */
  32 /*  teristic sets introduced by J. F. Ritt and developed by Wu Wen-  # */
  33 /*  tsun. The algorithms with variants implemented here  are on the  # */
  34 /*  basis of a generalization given by this author. Other modifica-  # */
  35 /*  tions are also made. For references, see                         # */
  36 /*  Ritt J. F. Differential Algebra. AMS, 1950.                      # */
  37 /*  Wang D. M. Characteristic Sets and Zero Structure of Polynomial  # */
  38 /*     Sets. Lecture Notes, RISC-LINZ, 1989.                         # */
  39 /*  Wu W. T. Basic Principles of Mechanical Theorem Proving in       # */
  40 /*     Elementary Geometries. J. Sys. Sci. & Math. Scis. 4 (1984)    # */
  41 /*     207-235; J. Automated Reasoning 2 (1986) 221-252.             # */
  42 /*===================================================================# */
  43 /* 1st change on May 29, 1991  */
  44 /* 2nd change on September 22, 1991 */
  45 /* 3rd change in December 1991 */
  46 /* 4th change in May 1992 */
  47 /* 5th change in October 1993 */
  48 /* 6th change in January 1994 */
  49 /* 7th change in April 1994 */
  50 /* 8th change in April 1995 */
  51 /* 9th change in August 1997 (for maxindset) */
  52 /* 10th change in February 1998 (for class and lvar) */
  53 /* Translated into macsyma by Dan Stanger */
  54 /* Info for maxima implementation: */
  55 /* Do not create variables with the name args, until the argument name */
  56 /* bug is fixed. */
  57
  58 load ("charsets_set.lisp");
  59 load ("charsets_flatten.lisp");
  60 load ("charsets_powers.lisp");
  61 load ("charsets_length.lisp");
  62
  63 /* load(grobner); */
  64
  65 define_variable(charsets_printlevel,0,any)$
  66 define_variable(charsets_qs,true,any,"Used as array")$
  67 define_variable(charsets_fact,1,any,"Stores a factor")$
  68 define_variable(charsets_ListOrSet,["[",set],any,"Used by charsets_operatorp")$
  69
  70 charsets_operatorp (e, l) :=
  71   if atom(e) then false else operatorp (e, l);
  72
  73
  74 define_variable(charsets_basset,[],any)$
  75 define_variable(charsets_brank,[],any)$
  76 define_variable(charsets_fwcharsetn,[],any)$
  77 define_variable(charsets_wcharsetn,[],any)$
  78 define_variable(charsets_con,[],any)$
  79 define_variable(charsets_fqcharsetn,[],any)$
  80 define_variable(charsets_wbasset,[],any)$
  81 define_variable(charsets_die,[],any)$
  82
  83
  84 /*#### Part 0. Definition of User Functions ##### */
  85 /* set of non-zero remainders of polys in ps wrt ascending set as */
  86 /*       user level function */
  87
  88 primpart(e,x):=part(content(e,x),2)$
  89
  90 /* This function is true if x is a name in the maple sense */
  91 namep(x):=is(symbolp(x) or subvarp(x))$
  92 /* This function is used to debug charsets_setify */
  93 setify1(x):=if listp(x) then charsets_setify(x) else charsets_setify([x])$
  94
  95 /* This is no more necessary, charsets_set.lisp has been fixed */
  96 /* charsets_union([x]):=listify(apply(union, maplist(lambda([y],setify(expand(y))),x)))$ */
  97
  98 /* This macro is used to deconstruct a list of 2 elements */
  99 ml2(u,v,l)::=buildq([u,v,l],
 100    (block([ml2_temp:l],u:first(ml2_temp),v:second(ml2_temp))))$
 101 ml3(u,v,w,l)::=buildq([u,v,w,l],(block([ml3_temp:l],
 102    u:first(ml3_temp),v:second(ml3_temp),w:third(ml3_temp))))$
 103 /* This function is used when there is a possibility that lcoeff could be
 104    called by -(x) */
 105 safeLcoeff(p,v):=
 106    if mapatom(p) then p /* what should lcoeff(x,y) be */
 107    else (
 108       if "-" = op(p) then -lcoeff(-p,v)
 109       else (
 110         if "*" = op(p) then safeLcoeff(first(p),v)*safeLcoeff(rest(p),v)
 111         else lcoeff(p,v)))$
 112
 113 /* this function is used to set up the charsets_qs variable */
 114 init_qs(ord):=block([l:length(ord)],
 115    if true = member(charsets_qs,arrays) then block([cl],
 116       cl:first(third(arrayinfo(charsets_qs))),
 117       if cl # l then array(charsets_qs,l))
 118    else array(charsets_qs,cl:l))$
 119
 120 /* this function is used to return the length of the charsets_qs variable */
 121 length_qs():= first(third(arrayinfo(charsets_qs)))$
 122 /* this function is used to flatten the charsets_qs array */
 123 flatten_qs(i,j):=block([o:[]],
 124    for ii:i thru j do (
 125       if listp(charsets_qs[ii]) then o:append(o,charsets_qs[ii])
 126       else o:endcons(charsets_qs[ii])),
 127    o
 128 )$
 129
 130 charsets_remset ( ps,as,ord):= block(
 131    [ _ind,i],
 132       if length ( ps)  < 1 or length ( as)  < 1 then (
 133             error ("no polynomials specified"))
 134          else if length ( ord)  < 1 then (
 135             error ("no indeterminates specified"))
 136          else if  not  listp(ord) then (
 137             error (ord, "must be a list")),
 138       if not every(map(namep, ord))
 139           then (
 140             error ("bad variable list")),
 141       _ind : 0,
 142       for i : 1 thru length ( as)  do
 143          if charsets_class ( as[i], ord)  <= _ind
 144              then (
 145                error ("second argument must be a non-contradictory (weak-, quasi-) ascending set"))
 146             else (
 147                _ind : charsets_class ( as[i], ord)),
 148       if  charsets_operatorp(ps, charsets_ListOrSet)
 149           then (
 150             if not every(map(lambda([p],polynomialp(p,ord,'ratnump)),ps))
 151                or
 152                not every(map(lambda([p],polynomialp(p,ord,'ratnump)),as))
 153                 then (
 154                   error ("input must be polynomials over Q in", ord)),
 155             charsets_remseta ( ps, as, ord))
 156          else (
 157             if not every(map(lambda([p],polynomialp(p,ord,'ratnump)),as))
 158                 then (
 159                   error ("input must be polynomials over Q in", ord)),
 160             charsets_premas ( ps, as, ord))
 161 )$
 162 /* the char set of poly set ps: user function */
 163
 164 charsets_charset ( ps,lst,[args]):= block( [ mset,ord,qs,y],
 165    if length ( ps)  < 1 then ( error ("no polynomials specified"))
 166    else if length ( lst)  < 1 then ( error ("no indeterminates specified"))
 167    else if 2 < length(args) then ( error ("too many arguments")),
 168    if not every(map(namep, lst)) then ( error ("bad variable list")),
 169    if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps)) then (
 170       error ("input must be polynomials over Q in", lst)),
 171    qs : charsets_setdifference( setify1(charsets_expand( ps)), setify1([ 0 ])),
 172    /* sort qs to try to get it in the same order as maple for testing */
 173    qs : sort(qs, charsets_setorderp),
 174    if listp(lst) then ( ord : lst)
 175    else ( ord : charsets_reorder ([ lst ], charsets_degord, qs)),
 176    if length(args) = 0 then (mset:charsetn, y:ord)
 177    else if length(args) = 1 then (mset:first(args), y:ord)
 178    else if length(args) = 2 then (mset:first(args), y:second(args)),
 179    if member(mset, setify1([charsetn, wcharsetn, qcharsetn, wbasset, qbasset,
 180          triset, trisetc, basset ])) then (
 181       charsets_charseta( qs, ord, concat(charsets_,mset)))
 182    else (
 183       error ("medial set must be one of `basset`,`wbasset`, `qbasset`, `charsetn`,`wcharsetn`,`qcharsetn`,`triset` and `trisetc`"))
 184 )$
 185
 186 /* modified from expand for lists */
 187 charsets_expand ( s) :=block(
 188       if listp(s)
 189           then (maplist(expand,s))
 190          else if charsets_setp(s) then (
 191              apply(setify1,maplist(expand,s)))
 192          else (
 193             expand ( s))
 194 )$
 195 /* the modified char set of poly set ps: user function */
 196
 197 charsets_mcharset ( ps,lst,[args]):= block( [ mset,ord,qs,y],
 198    if length ( ps)  < 1 then ( error ("no polynomials specified"))
 199    else if length ( lst)  < 1 then ( error ("no indeterminates specified"))
 200    else if 2 < length(args) then ( error ("too many arguments")),
 201    if not every(map(namep, lst)) then (error ("bad variable list")),
 202    if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps)) then (
 203       error ("input must be polynomials over Q in", lst)),
 204    qs : charsets_setdifference (  setify1(charsets_expand ( ps)),  setify1([ 0 ]) ),
 205    /* sort qs to try to get it in the same order as maple for testing */
 206    qs : sort(qs, charsets_setorderp),
 207    if listp(lst) then ( ord : lst)
 208    else ( ord : charsets_reorder (  (lst), charsets_degord, qs)),
 209    if length(args) = 0 then (mset:charsetn, y:ord)
 210    else if length(args) = 1 then (mset:first(args), y:ord)
 211    else if length(args) = 2 then (mset:first(args), y:second(args)),
 212    if member(mset,  setify1([ charsetn, wcharsetn, qcharsetn, wbasset,
 213            qbasset, triset, trisetc, basset ])) then (
 214       charsets_fcharseta ( qs, ord, concat(charsets_,mset)))
 215    else ( error ("medial set must be one of `basset`,`wbasset`,`qbasset`,`charsetn`,`wcharsetn`,`qcharsetn`,`triset` and `trisetc`"))
 216 )$
 217 /* the set of all nonconstant factors of initials of polys in as: user function */
 218
 219 charsets_iniset ( as,ord):=block(
 220    [ _ind,i],
 221    _ind : 0,
 222          if length ( as)  < 1 then (
 223             error ("no polynomials specified"))
 224          else if length ( ord)  < 1 then (
 225             error ("no indeterminates specified"))
 226          else if  not  listp(ord) then (
 227             error (ord, "must be a list")),
 228       if not every(map(namep, ord))
 229           then (
 230             error ("bad variable list")),
 231       if not every(map(lambda([p],polynomialp(p,ord,'ratnump)),as))
 232           then (
 233             error ("input must be polynomials over Q in", ord)),
 234       for i : 1 thru length ( as)  do
 235          if charsets_class ( as[i], ord)  <= _ind
 236              then (
 237                error ("first argument must be a non-contradictory (weak-, quasi-) ascending set"))
 238             else (
 239                _ind : charsets_class ( as[i], ord)),
 240       charsets_initialset ( charsets_expand ( as), ord)
 241 )$
 242 /* the char series of poly set ps: user function */
 243
 244 charsets_charser ( ps,lst,[args]):=block( [ mset,ord,qs, y],
 245    if length ( ps)  < 1 then ( error ("no polynomials specified"))
 246    else if length ( lst)  < 1 then ( error ("no indeterminates specified"))
 247    else if 2 < length(args) then ( error ("too many arguments")),
 248    if not every(map(namep, lst)) then ( error ("bad variable list")),
 249    if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps)) then (
 250       error ("input must be polynomials over Q in", lst)),
 251    qs : charsets_setdifference (  setify1(charsets_expand( ps)), setify1([ 0 ]) ),
 252    /* sort qs to try to get it in the same order as maple for testing */
 253    qs : sort(qs, charsets_setorderp),
 254    if listp(lst) then ( ord : lst)
 255    else ( ord : charsets_reorder (  [ (lst) ], charsets_degord, qs)),
 256    if length(args) = 0 then (mset:charsetn, y:ord)
 257    else if length(args) = 1 then (mset:first(args), y:ord)
 258    else if length(args) = 2 then (mset:first(args), y:second(args)),
 259    if member(mset,setify1([charsetn,wcharsetn,wbasset,trisetc,basset])) then (
 260       charsets_charseries ( qs, ord, concat(charsets_,mset)))
 261    else (
 262       error ("medial set must be one of 'basset','wbasset', 'charsetn','wcharsetn' and 'trisetc'"))
 263 )$
 264 /* the char series of poly set ps -- allowing to remove factors */
 265 /*       user function */
 266
 267 charsets_mcs ( ps,lst,[args]):=block(
 268    [ mset,ord,qs,y],
 269       if length ( ps)  < 1 then (
 270             error ("no polynomials specified"))
 271          else if length ( lst)  < 1 then (
 272             error ("no indeterminates specified"))
 273          else if 2 < length(args) then (
 274             error ("too many arguments")),
 275       if not every(map(namep, lst))
 276           then (
 277             error ("bad variable list")),
 278       if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps))
 279           then (
 280             error ("input must be polynomials over Q in", lst)),
 281       qs : charsets_setdifference (  setify1([ (charsets_expand ( ps)) ]),  setify1([ 0 ]) ),
 282       if listp(lst)
 283           then (
 284             ord : lst)
 285          else (
 286             ord : charsets_reorder (  [ (lst) ], charsets_degord, qs)),
 287       if length(args) = 0 then (mset:charsetn, y:ord)
 288       else if length(args) = 1 then (mset:first(args), y:ord)
 289       else if length(args) = 2 then (mset:first(args), y:second(args)),
 290       if member(mset,  setify1([ charsetn, wcharsetn, wbasset, trisetc, basset ]))
 291           then (
 292             charsets_fcharser ( qs, ord, concat(charsets_,mset)))
 293          else (
 294             error ("medial set must be one of 'basset','wbasset', 'charsetn','wcharsetn' and 'trisetc'"))
 295 )$
 296 /* the extended char series of poly set ps */
 297 /*      user function */
 298
 299 charsets_ecs ( ps,lst,[args]):=block( [ mset,ord,qs,y],
 300    if length ( ps)  < 1 then ( error ("no polynomials specified"))
 301    else if length ( lst)  < 1 then ( error ("no indeterminates specified"))
 302    else if 2 < length(args) then ( error ("too many arguments")),
 303    if not every(map(namep, lst)) then ( error ("bad variable list")),
 304    if listp(ps[1]) then (
 305       if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps[1])) then (
 306          error ("input must be polynomials over Q in", lst)))
 307       else if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps)) then (
 308          error ("input must be polynomials over Q in", lst)),
 309       if  charsets_operatorp(ps[1], charsets_ListOrSet) then (
 310          qs:[charsets_setdifference(setify1(charsets_expand(ps[1])),setify1([0])),ps[2]])
 311       else (
 312             qs : charsets_setdifference (setify1(charsets_expand ( ps)),setify1([ 0 ]))),
 313       if listp(lst) then ( ord : lst)
 314       else (
 315          if  charsets_operatorp(ps[1], charsets_ListOrSet) then (
 316             ord : charsets_reorder (  [ (lst) ], charsets_degord, qs[1]))
 317          else ( ord : charsets_reorder ( lst, charsets_degord, qs))),
 318       if length(args) = 0 then (mset:charsetn, y:ord)
 319       else if length(args) = 1 then (mset:first(args), y:ord)
 320       else if length(args) = 2 then (mset:first(args), y:second(args)),
 321       if member(mset,  setify1([ charsetn, wcharsetn, wbasset, trisetc, basset ]))
 322           then (
 323             charsets_excharser ( qs, ord, concat(charsets_,mset)))
 324          else (
 325             error ("medial set must be one of 'basset','wbasset', 'charsetn','wcharsetn' and 'trisetc'"))
 326 )$
 327 /* the extended char series of poly set ps -- allowing to remove factors */
 328 /*       user function */
 329
 330 charsets_mecs ( ps,lst,[args]):=block( [ mset,ord,qs,y],
 331    if length ( ps)  < 1 then ( error ("no polynomials specified"))
 332    else if length ( lst)  < 1 then ( error ("no indeterminates specified"))
 333    else if 2 < length(args) then ( error ("too many arguments")),
 334    if not every(map(namep, lst)) then ( error ("bad variable list")),
 335    if listp(ps[1]) then (
 336       if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps[1])) then (
 337          error ("input must be polynomials over Q in", lst)))
 338    else if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps)) then (
 339       error ("input must be polynomials over Q in", lst)),
 340    if  charsets_operatorp(ps[1], charsets_ListOrSet) then (
 341       qs:[sort(charsets_setdifference(setify1(charsets_expand (ps[1])),setify1([0])),
 342          charsets_setorderp),ps[2]])
 343    else (
 344       qs : sort(charsets_setdifference(setify1(charsets_expand(ps)), setify1([0])),
 345            charsets_setorderp)),
 346    if listp(lst) then ( ord : lst)
 347    else (
 348       if  charsets_operatorp(ps[1], charsets_ListOrSet) then (
 349          ord : charsets_reorder (  [ (lst) ], charsets_degord, qs[1]))
 350       else ( ord : charsets_reorder (  [ (lst) ], charsets_degord, qs))),
 351    if length(args) = 0 then (mset:charsetn, y:ord)
 352    else if length(args) = 1 then (mset:first(args), y:ord)
 353    else if length(args) = 2 then (mset:first(args), y:second(args)),
 354    if member(mset,  setify1([ charsetn, wcharsetn, wbasset, trisetc, basset ]))
 355           then ( charsets_fexcharser ( qs, ord, concat(charsets_,mset)))
 356    else (
 357       error ("medial set must be one of 'basset','wbasset', 'charsetn','wcharsetn' and 'trisetc'"))
 358 )$
 359 /* the irreducible char series of poly set ps: user function */
 360
 361 charsets_ics ( ps,lst,[args]):=block( [ mset,ord,qs,y],
 362    if length ( ps)  < 1 then ( error ("no polynomials specified"))
 363    else if length ( lst)  < 1 then ( error ("no indeterminates specified"))
 364    else if 2 < length(args) then ( error ("too many arguments")),
 365    if not every(map(namep, lst)) then ( error ("bad variable list")),
 366    if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps)) then (
 367        error ("input must be polynomials over Q in", lst)),
 368    qs : charsets_setdifference( setify1(charsets_expand ( ps)),  setify1([ 0 ]) ),
 369    if listp(lst) then ( ord : lst)
 370    else ( ord : charsets_reorder (  [ (lst) ], charsets_degord, qs)),
 371    if length(args) = 0 then (mset:charsetn, y:ord)
 372    else if length(args) = 1 then (mset:first(args), y:ord)
 373    else if length(args) = 2 then (mset:first(args), y:second(args)),
 374    if member(mset,  setify1([ charsetn, trisetc, basset ])) then (
 375        charsets_irrcharser ( qs, ord, concat(charsets_,mset)))
 376    else (error ("medial set must be one of 'basset','charsetn' and 'trisetc'"))
 377 )$
 378 /* the extended irreducible char series of poly set ps: user function */
 379
 380 charsets_eics ( ps,lst,[args]):=block( [ mset,ord,qs,y],
 381    if length ( ps)  < 1 then ( error ("no polynomials specified"))
 382    else if length ( lst)  < 1 then ( error ("no indeterminates specified"))
 383    else if 2 < length(args) then ( error ("too many arguments")),
 384    if not every(map(namep, lst)) then ( error ("bad variable list")),
 385    if  charsets_operatorp(ps[1], charsets_ListOrSet) then (
 386       if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps[1])) then (
 387          error ("input must be polynomials over Q in", lst)))
 388    else if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps)) then (
 389       error ("input must be polynomials over Q in", lst)),
 390    if  charsets_operatorp(ps[1], charsets_ListOrSet) then (
 391       qs:[sort(charsets_setdifference(setify1(charsets_expand (ps[1])),setify1([0])),
 392          charsets_setorderp),ps[2]])
 393    else (
 394       qs : sort(charsets_setdifference(setify1(charsets_expand(ps)), setify1([0])),
 395            charsets_setorderp)),
 396    if listp(lst) then ( ord : lst)
 397    else (
 398       if  charsets_operatorp(ps[1], charsets_ListOrSet) then (
 399          ord : charsets_reorder (  [ (lst) ], charsets_degord, qs[1]))
 400       else (
 401          ord : charsets_reorder (  [ (lst) ], charsets_degord, qs))),
 402    if length(args) = 0 then (mset:charsetn, y:ord)
 403    else if length(args) = 1 then (mset:first(args), y:ord)
 404    else if length(args) = 2 then (mset:first(args), y:second(args)),
 405    if member(mset,  setify1([ charsetn, trisetc, basset ])) then (
 406       charsets_exirrcharser(charsets_expand(qs),ord,concat(charsets_,mset)))
 407    else (
 408       error("medial set must be one of `basset`, `charsetn` and `trisetc`"))
 409 )$
 410 /* the quasi-irreducible char series of poly set ps: user function */
 411
 412 charsets_qics ( ps,lst,[args]):=block( [ mset,ord,qs,y],
 413    if length ( ps)  < 1 then ( error ("no polynomials specified"))
 414    else if length ( lst)  < 1 then ( error ("no indeterminates specified"))
 415    else if 2 < length(args) then ( error ("too many arguments")),
 416    if not every(map(namep, lst)) then ( error ("bad variable list")),
 417    if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps)) then (
 418       error ("input must be polynomials over Q in", lst)),
 419    qs : charsets_setdifference ( setify1(charsets_expand ( ps)),  setify1([ 0 ]) ),
 420    if listp(lst) then ( ord : lst)
 421    else ( ord : charsets_reorder (  [ (lst) ], charsets_degord, qs)),
 422    if length(args) = 0 then (mset:charsetn, y:ord)
 423    else if length(args) = 1 then (mset:first(args), y:ord)
 424    else if length(args) = 2 then (mset:first(args), y:second(args)),
 425    if member(mset,setify1([charsetn,wcharsetn,wbasset,trisetc,basset])) then (
 426       charsets_qirrcharser ( qs, ord, concat(charsets_,mset)))
 427    else (
 428       error ("medial set must be one of `basset`,`wbasset`, `charsetn`,`wcharsetn` and `trisetc`"))
 429 )$
 430 /* factorize poly f over algebraic number field with minimal polys in as */
 431 /*       wrt ord: user function */
 432 /*       used for linear transformation */
 433
 434 charsets_cfactor ( f,[args]):=block( [ _ind,inda,ff,i,as,ord,last],
 435    if length(args) = 0 then ( return ( factor ( f) )),
 436    if length(args) # 2 then ( error ("improper number of arguments"))
 437    else (as:first(args),ord:second(args)),
 438    if length ( as)  < 1 then ( error ("no polynomials specified"))
 439    else if length ( ord)  < 1 then ( error ("no indeterminates specified"))
 440    else if  not  listp(ord) then ( error (ord, "must be a list"))
 441    else if not every(map(namep, ord)) then ( error ("bad variable list"))
 442    else if not every(map(lambda([p],polynomialp(p,ord,'ratnump)),as)) then (
 443       error ("input must be polynomials over Q in", ord)),
 444    ff : num( f),
 445    _ind : 0,
 446    for i : 1 thru length ( as)  do (
 447       inda : charsets_class ( as[i], ord),
 448       if inda <= _ind then (
 449          error ("second argument must be a non-contradictory ascending set"))
 450       else ( _ind : inda)),
 451    if 1 < charsets_printlevel then (
 452       print("Warning: be sure the ascending set is irreducible")),
 453    if charsets_class( ff, ord) <= charsets_class( as[length( as)], ord) then (
 454       factor ( f))
 455    else (
 456       last:apply("+", maplist( lambda([p],charsets_degreel( p, ord)), as)),
 457       charsets_das :if charsets_degreel ( ff, ord)  < last then
 458          [ -1, 1, -2, 2, -3, false] else [ 1, -1, 2, -2, -3, false],
 459       last:factor ( f),
 460       if "-" = op(last) then last:-last,
 461       charsets_cfactorsub ( last, as, ord))
 462 )$
 463 /* prepare a list of triangular forms from poly set ps: user function */
 464 /* I am not sure if I got the argument handling correct */
 465 /* In fact the argument handling is incorrect. Either it is triser(ps,ord) */
 466 /* where ord is a list, or triser(ps,lst,ord) where lst is the unordered set */
 467 /* of variables and ord the order computed by charsets_reorder (lst, charsets_degord, qs) */
 468 /* and automatically filled */
 469
 470 charsets_triser ( ps,[args]):=block( [ i,ord,qs, lst,y],
 471    if length ( ps)  < 1 then ( error ("no polynomials specified"))
 472    else if 2 < length(args) then ( error ("too many arguments"))
 473    else if length(args) >= 1 then (
 474       lst:first(args),
 475       if length ( lst) < 1 then ( error ("no indeterminates specified")),
 476       if not every(map(namep, lst)) then ( error ("bad variable list"))),
 477    if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps)) then (
 478       error ("input must be polynomials over Q in", lst)),
 479    qs : charsets_setdifference( setify1(charsets_expand ( ps)), setify1([ 0 ])),
 480    if length(args) < 1 then (
 481       ord : charsets_reorder ( maplist(listofvars, ps), charsets_degord, qs))
 482    else if listp(lst) then ( ord : lst)
 483    else ( /* this is never executed? as length(lst) will fail above */
 484       ord : charsets_reorder ( lst, charsets_degord, qs)),
 485    if length(args) = 1 then y : ord else y:second(args),
 486    charsets_trisersub ( qs, ord)
 487 )$
 488 /* solve a set of poly eqs ps=0: user function */
 489 /* Arg handling is copied from triser */
 490
 491 charsets_csolve ( ps,[args]):=block( [ i,ord,qsi,qs, lst,y],
 492    if length ( ps)  < 1 then ( error ("no polynomials specified"))
 493    else if 2 < length(args) then ( error ("too many arguments"))
 494    else if length(args) >= 1 then (
 495       lst:first(args),
 496       if length ( lst) < 1 then ( error ("no indeterminates specified")),
 497       if not every(map(namep, lst)) then ( error ("bad variable list"))),
 498    if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps)) then (
 499       error ("input must be polynomials over Q in", lst)),
 500    qs : charsets_setdifference( setify1(charsets_expand ( ps)), setify1([ 0 ])),
 501    if length(args) < 1 then (
 502       ord : charsets_reorder ( maplist(listofvars, ps), charsets_degord, qs))
 503    /* else if listp(lst) then ( ord : lst) Needed if sets are not lists */
 504    else ( ord : charsets_reorder( lst, charsets_degord, qs)),
 505    if length(args) = 1 then y : ord else y:second(args),
 506    /* sort to try to get things in maple order */
 507    qs : sort (qs, charsets_lenord),
 508    qsi :  setify1( charsets_trisersub ( qs, ord)),
 509    if qsi =  setify1(setify1([])) then ( setify1([]))
 510    else setify1( maplist(lambda([q],charsets_solveas( q, ord)), qsi))
 511 )$
 512 /* irreducible decomposition of algebraic variety defined by ps */
 513 /*      user function */
 514
 515 charsets_ivd ( ps,lst,[args]):=block( [ mset,ord,qs,y],
 516    if length ( ps)  < 1 then ( error ("no polynomials specified"))
 517    else if length ( lst)  < 1 then ( error ("no indeterminates specified"))
 518    else if 2 < length(args) then ( error ("too many arguments")),
 519    if not every(map(namep, lst)) then ( error ("bad variable list")),
 520    if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps)) then (
 521       error ("input must be polynomials over Q in", lst)),
 522    qs : charsets_setdifference( setify1(charsets_expand ( ps)), setify1([ 0 ]) ),
 523    if listp(lst) then ( ord : lst)
 524    else ( ord : charsets_reorder (  [ (lst) ], charsets_degord, qs)),
 525    if length(args) = 0 then (mset:charsetn, y:ord)
 526    else if length(args) = 1 then (mset:first(args), y:ord)
 527    else if length(args) = 2 then (mset:first(args), y:second(args)),
 528    if member(mset,  setify1([ charsetn, trisetc, basset ])) then (
 529       charsets_irrvardec ( qs, ord, concat(charsets_,mset)))
 530    else (
 531       error ("medial set must be one of `basset`,`charsetn` and `trisetc`"))
 532 )$
 533 /* primary decomposition of polynomial ideal generated by ps */
 534 /*      user function */
 535
 536 charsets_pid ( ps,lst,[args]):=block( [ mset,ord,qs,y],
 537    if length ( ps)  < 1 then ( error ("no polynomials specified"))
 538    else if length ( lst)  < 1 then ( error ("no indeterminates specified"))
 539    else if 2 < length(args) then ( error ("too many arguments")),
 540    if not every(map(namep, lst)) then ( error ("bad variable list")),
 541    if not every(map(lambda([p],polynomialp(p,lst,'ratnump)),ps)) then (
 542       error ("input must be polynomials over Q in", lst)),
 543    qs : charsets_setdifference (  setify1(charsets_expand ( ps)),  setify1([ 0 ]) ),
 544    if listp(lst) then ( ord : lst)
 545    else ( ord : charsets_reorder (  [ (lst) ], charsets_degord, qs)),
 546    if length(args) = 0 then (mset:charsetn, y:ord)
 547    else if length(args) = 1 then (mset:first(args), y:ord)
 548    else if length(args) = 2 then (mset:first(args), y:second(args)),
 549    if member(mset,  setify1([ charsetn, trisetc, basset ])) then (
 550       charsets_pridealdec ( qs, ord, concat(charsets_,mset)))
 551    else (
 552       error ("medial set must be one of `basset`, `charsetn` and `trisetc`"))
 553 )$
 554 /*#### Part I. Routines for Computing Characteristic Sets ##### */
 555 /* the class of poly f wrt variable ordering ord */
 556
 557 charsets_class ( f,ord):=block(
 558    [ ii:0],
 559       for i : length ( ord)  thru 1 step  -  1 do
 560          if member(ord[i], listofvars( f))
 561              then (
 562                return ( ii:i)),
 563       ii
 564 )$
 565 /* the leading variable of poly f wrt variable ordering ord    */
 566
 567 charsets_lvar ( f,ord):=block([ i,r:0],
 568       for i : length ( ord)  thru 1 step  -  1 do
 569          if member(ord[i], listofvars ( f)) then (
 570                return ( r:ord[i])),
 571       if r # 0 then return(r),
 572       if 1 < charsets_printlevel then (
 573          print("Warning: lvar is called with constant")),
 574       0
 575 )$
 576 /* the index set of a poly (or a poly set f) wrt ord */
 577
 578 charsets_index( f,ord):=block( [ i,g,ng,O],
 579    if listp(f) then (
 580        map(lambda([p],charsets_index( p, ord)), f))
 581    else if /* WARNING: listp should be charsets_setp */listp(f) then (
 582        apply(setify1,map(lambda([p],charsets_index( p, ord)), f)))
 583    else (
 584       g : expand ( f),
 585       ng : charsets_terms ( g),
 586       if charsets_class ( g, ord)  = 0 then (
 587          i : part(  [ (listofvars ( g)), O ], 1))
 588       else (
 589          i : charsets_lvar ( g, ord)),
 590       g : degree ( g, i),
 591       [ng, i, g])
 592 )$
 593 /* number of terms in g */
 594
 595 charsets_terms ( g):=block(
 596       if "+" = inpart(g, 0)
 597           then ( length ( g))
 598          else ( 1)
 599 )$
 600 /* the initial of poly p wrt ord */
 601
 602 charsets_initial ( p,ord):=block(
 603    [ f,last],
 604       f : expand ( p),
 605       if charsets_class ( f, ord)  = 0 then 1
 606          else (
 607             last:lcoeff ( f, charsets_lvar ( f, ord) ),
 608             /* rncombine worked here but is not implemented in maxima */
 609             num(xthru(multthru( last / lcoeff ( last,ord)))))
 610 )$
 611 /* modified rank of two polys: comparing further the rank */
 612 /*     of initials when f and g have same rank */
 613
 614 charsets_mrank ( f,g,ord):=block(
 615    [ cf,cg],
 616       cf : charsets_class ( f, ord),
 617       cg : charsets_class ( g, ord),
 618       if cf = 0 then ( true)
 619          else if cf < cg then ( true)
 620          else if cf = cg then (
 621             cf : charsets_degreel ( f, ord),
 622             cg : charsets_degreel ( g, ord),
 623             if cf < cg then ( true)
 624                else if cf = cg then (
 625                   charsets_mrank ( charsets_initial ( f, ord), charsets_initial ( g, ord), ord))
 626                else ( false))
 627          else ( false)
 628 )$
 629 /* modified rank of two polys: comparing further the rank of */
 630 /*     initials, the terms of initials and the terms of f and g */
 631 /*     when they are the same */
 632
 633 charsets_rank ( f,g,ord):=block(
 634    [ _ind,find,cf,cg],
 635       ml2(find,_ind, charsets_subrank ( f, g, ord)),
 636       if find and _ind = 1 then (
 637             cf : charsets_polylength ( expand ( charsets_initial ( f, ord) ) ),
 638             cg : charsets_polylength ( expand ( charsets_initial ( g, ord) ) ),
 639             if cf < cg
 640                 then ( true)
 641                else if cf = cg then (
 642                   if charsets_polylength ( expand ( f) )  < charsets_polylength ( expand ( g) )
 643                       then ( true)
 644                      else ( false))
 645                else ( false))
 646          else ( find)
 647 )$
 648 /* subroutine for rank */
 649
 650 charsets_subrank ( f,g,ord):=block(
 651    [ cf,cg,_ind:0],
 652       cf : charsets_class ( f, ord),
 653       cg : charsets_class ( g, ord),
 654       if cf = 0 then (
 655             if cg = 0 then (_ind : 1), [true,_ind])
 656          else if cf < cg then ( [true,_ind])
 657          else if cf = cg then (
 658             cf : charsets_degreel ( f, ord),
 659             cg : charsets_degreel ( g, ord),
 660             if cf < cg then ( [true,_ind])
 661                else if cf = cg then (
 662                   charsets_subrank ( charsets_initial ( f, ord),
 663                        charsets_initial ( g, ord), ord))
 664                else ( [false,_ind]))
 665          else ( [false,_ind])
 666 )$
 667 /* the rank of two polys with same classes:  */
 668 /*        used for computing tiangular form  */
 669
 670 charsets_trank ( f,g,ord):=block(
 671    [ cf,cg],
 672       cf : charsets_degreel ( f, ord),
 673       cg : charsets_degreel ( g, ord),
 674       if cf < cg
 675           then (
 676             true)
 677
 678          else if cf = cg then (
 679             charsets_mrank ( charsets_initial ( f, ord), charsets_initial ( g, ord), ord))
 680          else (
 681             false)
 682 )$
 683 /* modified pseudo division: I1^s1...Ir^sr*uu = q*vv + r, */
 684 /*    where I1, ..., I_r are all distinct factors of lcoeff(vv,x) */
 685 /*    and s1, ..., sr are chosen to be the smallest  */
 686 charsets_prem( uu,vv,x):=block(
 687    if charsets_operatorp(vv,charsets_ListOrSet) = true then (
 688       if length(vv) > 1 then error("charsets_prem length vv > 1, not implemented", vv)
 689       else vv:first(vv)),
 690    if charsets_operatorp(uu,charsets_ListOrSet) = true then
 691       map(lambda([u],charsets_prem1(u,vv,x)),uu)
 692    else
 693       charsets_prem1(uu,vv,x)
 694 )$
 695 charsets_prem1( uu,vv,x):=block( [ r,v,dr,dv,l,t,lu,lv,g,di],
 696    if integerp(vv / x) then ( subst(x = 0,uu))
 697    else (
 698       r : expand ( uu),
 699       dr : degree ( r, x),
 700       v : expand ( vv),
 701       dv : degree ( v, x),
 702       if dv <= dr then (
 703                   l : coeff(v,x,dv),
 704                   v : expand ( v - l * x ^ dv))
 705       else ( l : 1),
 706       if 1 < charsets_printlevel and 500 < length ( r)  + length ( v) then (
 707          print("pseudo-division:",[length(r),x,dr], [ length ( v), x, dv ])),
 708       while (dv <= dr and r # 0) do (
 709          g:gcd ( lu:l, lv:coeff(r,x,dr)),
 710          lu:first(di:divide(lu,g)),
 711          if second(di) # 0 then error("prem1 divide by gcd yielded remainder!",
 712              di,lu,lv),
 713          lv:first(di:divide(lv,g)),
 714          if second(di) # 0 then error("prem1 divide by gcd yielded remainder!",
 715              di,lu,lv),
 716          t : expand ( x ^ ( ( dr - dv)) * v * lv),
 717          if dr = 0 then ( r : 0)
 718          else ( r : subst(x ^ dr = 0, r)),
 719          r : expand ( lu * r)  - t,
 720          dr : degree ( r, x)),
 721    r)
 722 )$
 723 /* pseudo remainder of poly f wrt ascending set as */
 724
 725 charsets_premas ( f,as,ord):=block([ remd,i],
 726    remd : f,
 727    for i : length ( as)  thru 1 step  -  1 do
 728       remd : charsets_prem ( remd, as[i], charsets_lvar ( as[i], ord) ),
 729    if remd # 0 then ( num(xthru(multthru( remd / lcoeff( remd,/* reverse */(ord))))))
 730    else ( 0)
 731 )$
 732 /* set of non-zero remainders of polys in ps wrt ascending set as */
 733
 734 charsets_remseta ( ps,as,ord):=block([retval],
 735    retval:sort(charsets_setdifference (
 736       setify1(maplist(lambda([p],charsets_premas( p, as, ord)),ps)),
 737       setify1([ 0 ])),charsets_setorder_reversep),
 738    retval
 739 )$
 740 /* pseudo remainder of poly f wrt ascending set as -- version b */
 741
 742 charsets_premasb ( f,as,ord):=block(
 743    [ remd,i,l],
 744    remd : f,
 745    if 1 < length ( as)  then (
 746       for i : length ( as)  thru 2 step  -  1 do
 747          remd : charsets_prem ( remd, as[i], charsets_lvar ( as[i], ord))),
 748    l:divide ( remd, primpart(as[1],charsets_lvar(as[1],ord))),
 749    if first(l)=0  then ( remd : 0)
 750    else ( remd : charsets_prem ( remd, as[1], charsets_lvar ( as[1], ord))),
 751    if remd # 0 then ( num(xthru(multthru( remd / lcoeff( remd,/* reverse */(ord))))))
 752    else (0)
 753 )$
 754 /* set of non-zero remainders of polys in ps wrt ascending set as -- version b */
 755
 756 charsets_remsetb ( ps,as,ord):=block([last],
 757    last:maplist(lambda([p],charsets_premasb( p, as, ord)),ps),
 758    charsets_setdifference( setify1(if listp(last) then last else [last]),setify1([0]))
 759 )$
 760 /* reorder the list ord of variables wrt poly set ps */
 761
 762 charsets_reorder ( ord,p,ps):=block([last],
 763    last:charsets_reordera ( ord, ps),
 764    append(charsets_reorderb(charsets_setdifference(setify1(ord),setify1(last)),p,ps),
 765       last)
 766 )$
 767 /* subroutine for reorder: first criterion */
 768
 769 charsets_reordera ( ord,ps):=block([ qs,pp,orb,i,last],
 770    if length ( ps)  = 0 then ( ord)
 771    else (
 772       last:[],
 773       qs :  setify1([ (ps) ]),
 774       orb :  setify1([ (ord) ]),
 775       for i in orb do (
 776          pp : charsets_deg0 ( ps, i),
 777          if length ( pp)  = 1 then (
 778             return (last:cons(i,charsets_reordera(charsets_setdifference(orb,setify1(i)),
 779                charsets_setdifference( qs, pp) ))))),
 780    last)
 781 )$
 782
 783 /* subroutine for reorder -- modified from sort: second criterion */
 784 /* local works in maxima, but not in macsyma */
 785 charsets_reorderb ( l,p,ps):=block([ n,tn,gap,ggap,i,j,temp,v],
 786    n : length ( l),
 787    local(v),array(v,n),
 788    tn : p,
 789    for i : 1 thru n do v[i - 1] : l[i],
 790    for ggap:(gap:4) while ggap <= n do (gap:(ggap : 3*ggap+1)),
 791    gap : quotient( gap, 3),
 792    while (0 < gap) do (
 793       for i : gap thru n - 1 do (
 794          temp : v[i],
 795          j:i - gap,
 796          while (j >= 0) do (
 797             if tn ( v[j], temp, ps) then return(true),
 798             v[j + gap] : v[j],
 799             j:j-gap
 800          ),
 801          v[j + gap] : temp
 802       ),
 803       gap : quotient( gap, 3)
 804    ),
 805    temp:makelist(v[i],i,0,n-1),
 806    temp
 807 )$
 808 /* determine the order between x and y wrt ps */
 809
 810 charsets_degord ( x,y,ps):=block(
 811    if part(charsets_degpsmax(ps,y),2) < part(charsets_degpsmax(ps,x),2)
 812           then true
 813    else if part(charsets_degpsmax(ps,x),2) < part(charsets_degpsmax(ps,y),2)
 814           then false
 815    else if part(charsets_degpsmax(ps,y),1) < part(charsets_degpsmax(ps,x),1)
 816           then true
 817    else if part(charsets_degpsmax(ps,x),1) < part(charsets_degpsmax(ps,y),1)
 818           then false
 819    else if part(charsets_degpsmin(ps,x),2) < part(charsets_degpsmin(ps,y),2)
 820           then true
 821    else if part(charsets_degpsmin(ps,y),2) < part(charsets_degpsmin(ps,x),2)
 822           then false
 823    else if part(charsets_degpsmin(ps,y),1) < part(charsets_degpsmin(ps,x),1)
 824           then true
 825    else if part(charsets_degpsmin(ps,x),1) < part(charsets_degpsmin(ps,y),1)
 826           then false
 827    else if part(charsets_deg1(ps,y),1) < part(charsets_deg1(ps,x),1) then true
 828    else if part(charsets_deg1(ps,x),1) < part(charsets_deg1(ps,y),1) then false
 829    else if part(charsets_deg1(ps,y),2) < part(charsets_deg1(ps,x),2) then true
 830    else false
 831 )$
 832 /* the maximal degree of polys in qs wrt x  */
 833 /*      and the number of polys having this degree  */
 834
 835 charsets_degpsmax ( qs,x):=block(
 836    [ i,m,mm,ps],
 837       /* Is this charsets_expand? */
 838       ps : expand ( qs),
 839       m : apply(max,maplist(lambda([p],degree(p, x)),ps)),
 840       mm : 0,
 841       for i : 1 thru length ( ps)  do
 842          if degree ( ps[i], x)  = m then ( mm : mm + m),
 843        [ mm, m ]
 844 )$
 845 /* the minimal non-zero degree of polys in qs wrt x  */
 846 /*      and the number of polys having this degree  */
 847
 848 charsets_degpsmin ( qs,x):=block( [ i,m,mm,ps,last],
 849    ps : expand ( qs),
 850    last:charsets_setdifference( setify1(map(lambda([p],degree( p, x)),ps)),
 851           setify1([ 0 ]) ),
 852    if last =  setify1([]) then ( m : 0)
 853    else m : apply(min, last),
 854    mm : 0,
 855    for i : 1 thru length ( ps)  do
 856       if degree ( ps[i], x)  = m then ( mm : mm + m),
 857    [ mm, m ]
 858 )$
 859 /* determine if ps has one and only one poly involving x */
 860
 861 charsets_deg0 ( ps,x):=block( [ i,ms],
 862    ms :  setify1([]),
 863    for i in ps do (
 864       if length( ms)  >= 2 then return(false),
 865       if member(x, listofvars( i)) then ms : charsets_adjoin2( i, ms)),
 866    ms
 867 )$
 868 /* the minimal total degree of lcoeffs of polys in PS wrt x  */
 869 /*      and the minimal number of terms of those lcoeffs */
 870 /*      I think this function was memoized */
 871
 872 charsets_deg1 ( _PS,x):=block( [ i,ps,qs,k],
 873    ps : expand ( _PS),
 874    qs :  setify1([]),
 875    k : part( charsets_degpsmin( ps, x), 2),
 876    for i : 1 thru length ( ps)  do (
 877       if degree( ps[i], x) = k then qs : charsets_adjoin2(lcoeff( ps[i], x),qs)),
 878    [ apply(min, maplist(lambda([u],degree( u, listofvars( u))),qs)),
 879      apply(min, maplist(lambda([u],charsets_polylength( expand (u))),qs))]
 880 )$
 881 /* search an element with lowest rank in ps */
 882 /*      and assign the rest of polys to qs */
 883
 884 charsets_sort ( ps,theRankFunction,ord):=block(
 885    [ l,i,qs],
 886       if length ( ps)  = 1 then ( [ps[1],[]])
 887          else (
 888             l : ps[1],
 889             qs :  [],
 890             for i : 2 thru length ( ps)  do
 891                if theRankFunction( ps[i], l, ord)
 892                    then (
 893                      qs :  cons( l, qs),
 894                      l : ps[i])
 895                   else (
 896                      qs : cons( ps[i], qs)),
 897             [l,qs])
 898 )$
 899 /* the difference of two lists */
 900
 901 charsets_minus ( ps,qs):=block(
 902    charsets_setdifference (  setify1(ps),  setify1(qs))
 903 )$
 904 /* the product of all elements in a list */
 905
 906 charsets_prod ( ps):=block( [ i],
 907    product ( ps[i], i, 1, length ( ps) )
 908 )$
 909 charsets_degree ( f,x):=block([c:facsum( f, x)],
 910    if mapatom(x) then charsets_degree(c,[x]),
 911    if listp(f) then lmax(maplist(lambda([u],degree(u,x)),c))
 912    else charsets_degree([f],x)
 913 )$
 914 charsets_degreel ( f,ord):=block([last],
 915       last:expand ( f),
 916       degree ( last, charsets_lvar ( last, ord) )
 917 )$
 918 /* the macro which creates the set functions */
 919 makeBasicSetFunction(theName,theRank)::=buildq([theName,theRank],
 920    theName( ps,ord):=block([ qs,qs1,i,b],
 921       if length ( ps)  < 2 then ( ps)
 922       else (
 923          ml2(b,qs1, charsets_sort ( ps, charsets_rank, ord)),
 924          qs : [],
 925          if 0 < charsets_class ( b, ord) then (
 926              for i in qs1 do (
 927                 if theRank( i, b, ord) then (
 928                     qs : cons( i, qs))))
 929          else (
 930             return ([ b ])),
 931          cons( b, theName( qs, ord)))
 932   ))$
 933 /* the basic set of poly set ps */
 934 makeBasicSetFunction(charsets_basset,charsets_brank)$
 935 /* the weak-basic set of poly set ps */
 936 makeBasicSetFunction(charsets_wbasset,charsets_wbrank)$
 937 /* the quasi-basic set of poly set ps */
 938 makeBasicSetFunction(charsets_qbasset,charsets_qbrank)$
 939 /* subroutines for basset, wbasset and qbasset */
 940
 941 charsets_brank ( a,b,ord):=block(
 942    if charsets_degree( a, charsets_lvar( b, ord)) < charsets_degreel( b, ord)
 943           then ( true)
 944    else ( false)
 945 )$
 946
 947 charsets_wbrank ( a,b,ord):=block(
 948    if charsets_class( b, ord) < charsets_class( a, ord) and
 949       charsets_brank ( charsets_initial ( a, ord), b, ord) then ( true)
 950    else ( false)
 951 )$
 952
 953 charsets_qbrank ( a,b,ord):=block(
 954    if charsets_class( b, ord) < charsets_class( a, ord) then true else false
 955 )$
 956
 957 /* the char set of poly set ps */
 958 /* add univariate gcd on May 5, 2000 */
 959 charsets_charseta ( ps,ord,medset):=block(
 960    [ cs,rs,l,med],
 961    if length ( ps)  < 2 then (ps)
 962    else (
 963       ps:charsets_unigcd(ps,ord),
 964       if medset = charsets_qcharsetn then (
 965          charsets_with :  setify1([]),
 966          cs : charsets_qcharsetn ( ps, ord, charsets_remsetaA))
 967       else ( cs : medset ( ps, ord)),
 968       if 0 < charsets_class ( cs[1], ord) then (
 969          if member(medset,setify1([charsets_wbasset,charsets_qbasset,
 970                 charsets_basset])) then (
 971             rs:charsets_remseta(charsets_setdifference(setify1(ps),setify1(cs)),cs,ord))
 972          else if medset = charsets_qcharsetn and charsets_checkwith(
 973                 charsets_with, charsets_initialset1 ( cs, ord) )  then (
 974             if 1 < charsets_printlevel then (
 975                print("Strategy for 0 remainder succeed")),
 976                return ( cs))
 977             else (
 978               if medset = charsets_qcharsetn and 1 < charsets_printlevel then (
 979                 print("Strategy for 0 remainder failed")),
 980               block([temp:charsets_canonlt,sd],
 981                 charsets_canonlt:charsets_setorder_reversep,
 982                 sd:charsets_setdifference(
 983                   setify1(if listp(ps) then ps else [ps]),
 984                   setify1(if listp(cs) then cs else [cs])),
 985                 sd:sort(sd,lambda([x,y],not charsets_lenord(x,y))),
 986                 rs : charsets_remsetb( sd, cs, ord),
 987                 charsets_canonlt:temp
 988                 )))
 989           else (
 990          return (  [ 1 ] )),
 991        if rs =  setify1([]) then (
 992          changesign(maplist(lambda([c],num(c/lcoeff(expand(c),/* reverse */(ord)))),cs)))
 993        else charsets_charseta( charsets_union( rs, cs, ps), ord, medset))
 994 )$
 995 /*
 996 changesign(l):=block([i],
 997    i:read("changesign 0 or item #"),
 998    if i # 0 then l[i]:-l[i],
 999    l)$
1000 */
1001 changesign(l):=l$
1002 charsets_unigcd(PS,ord):=block([_ps,p,c],
1003      _ps :charsets_setify([]),
1004      for p in PS do (
1005          c : charsets_class(p,ord),
1006          if c = 0 then return(charsets_setify([1]))
1007          else (if c = 1 then _ps : charsets_union(_ps, charsets_setify( [p])))
1008     ),
1009      if length(_ps) < 2 then PS
1010      else charsets_union(charsets_setdifference(PS,_ps), setify1(apply(charsets_GCD,_ps)))
1011 )$
1012
1013 charsets_GCD([args]) :=
1014      if length(args) = 2 then apply('gcd,args)
1015      else gcd(first(args),apply(charsets_GCD(rest(args))))$
1016
1017 /* the modified char set of poly set ps */
1018
1019 charsets_fcharseta ( ps,ord,medset):=block(
1020    [ csf,fset],
1021       ml2(csf,fset, charsets_fcharsetsub ( charsets_nopower ( ps), ord, medset,
1022           [  setify1([]), listofvars ( ps)  ])),
1023    [csf,factors_removed:fset[1]]
1024
1025 )$
1026 /* the main subroutine for fcharseta */
1027 /* add univariate gcd and change removecont on May 5, 2000 */
1028 define_variable(charsets_fcharsetsub_alist,
1029    maplist(lambda([x], concat(charsets_,x)=concat(charsets_f,x)),
1030       [charsetn, wcharsetn, qcharsetn, triset, trisetc]),any_check,
1031       "Used to convert medset to fmedset, BTW any_check didnt work")$
1032
1033 charsets_fcharsetsub ( _ps,ord,medset,fset1):=block(
1034    [ ps,cs,rs,l,fset3,fset2,ts,fmedset,med,fset,last],
1035    if length ( _ps)  < 2 then ( [_ps,fset1])
1036    else (
1037       ps : charsets_unigcd(_ps,ord),
1038       fmedset:assoc(medset,charsets_fcharsetsub_alist),
1039       if fmedset # false then (
1040          if fmedset = charsets_fqcharsetn then (
1041             charsets_with :  setify1([]),
1042             ml2(cs,fset2,charsets_fqcharsetn(ps, ord, fset1, charsets_remsetaA))
1043         )
1044          else (ml2(cs,fset2, fmedset ( ps, ord, fset1))),
1045          last:if length( cs) > 1 then cs[1] else [],
1046          if 0 < charsets_class ( last, ord) then (
1047             if fmedset = charsets_fqcharsetn and
1048                charsets_checkwith ( charsets_with,
1049                     charsets_union ( charsets_initialset1 ( cs, ord), fset2[1])) then (
1050                if 1 < charsets_printlevel then (
1051                   print("Strategy for 0 remainder succeed")),
1052                return ( [cs,fset2]))
1053             else (
1054                if medset = charsets_qcharsetn and 1 < charsets_printlevel then (
1055                   print("Strategy for 0 remainder failed")),
1056                   rs : charsets_remsetb(charsets_setdifference(setify1(ps),setify1(cs)),
1057                      cs, ord)))
1058          else (
1059             return ( [[ 1 ],fset2] )),
1060          if rs =  setify1([]) then (
1061             [map(lambda([c],num(c/lcoeff(expand(c),/* reverse */(ord)))),cs),fset2])
1062          else (
1063             charsets_fcharsetsub(charsets_union(rs,cs, ps),ord,medset,fset2)))
1064       else (
1065          if medset = [] then cs: ps else  cs : medset ( ps, ord),
1066          fset2 : [ fset1[1], charsets_union( fset1[2], charsets_initialset1(cs, ord))],
1067          if 0 < charsets_class ( cs[1], ord) then (
1068             rs:charsets_remseta(charsets_setdifference(setify1(ps),setify1(cs)),cs,ord),
1069             if length(rs) # 0  then (
1070                ml2(rs,ts,charsets_removecont(rs,ord)),
1071                ml2(rs,fset3,charsets_removefactor ( rs, ord, fset2)),
1072                fset3:[charsets_union(fset3[1],ts),fset3[2]])
1073             else ( rs:[], fset3:fset2))
1074          else (
1075             return ( [[ 1 ],fset2])),
1076          if rs =  []  then (
1077             [map(lambda([c],num(c/lcoeff(expand(c),/* reverse */(ord)))),cs),fset3])
1078          else
1079             charsets_fcharsetsub ( charsets_union ( rs, cs, ps), ord,
1080                    medset, fset3)))
1081 )$
1082 /*## The following few routines implement a strategy for speeding-up  */
1083 /*## the computation of charsets by remembering all appearing initials */
1084 /*## in the quasi-sense. */
1085
1086 charsets_premA ( uu,vv,x):=block( [ r,v,dr,dv,l,t,lu,lv,g],
1087    if integerp(vv / x) then ( subst( x = 0,uu))
1088    else (
1089       r : expand ( uu),
1090       dr : degree ( r, x),
1091       v : expand ( vv),
1092       dv : degree ( v, x),
1093       if dv <= dr then (
1094          l : coeff(v,x,dv),
1095          v : expand ( v - l * x ^ dv))
1096       else ( l : 1),
1097       if 1 < charsets_printlevel and 500 < length ( r)  + length ( v) then (
1098          print("pseudo-division:, [ length ( r), x, dr ] ",[length(v),x,dv])),
1099       while (dv <= dr and r # 0) do (
1100          g:gcd ( lu:l, lv:coeff(r,x,dr)), lu:lu/g, lv:lv/g,
1101          t : expand ( x ^ ( ( dr - dv)) * v * lv),
1102          if dr = 0 then ( r : 0)
1103          else ( r : subst(x ^ dr = 0, r)),
1104          r : expand ( lu * r)  - t,
1105          if not ratnump(lu) and /* WARNING: listp should be charsets_setp */listp(charsets_with)
1106                    then (
1107             charsets_with : charsets_union ( charsets_with,  setify1([ lu ]) )),
1108             dr : degree ( r, x)
1109       ),
1110       r)
1111 )$
1112 charsets_remsetaA ( ps,as,ord):=block( [ i],
1113    charsets_setdifference(setify1([map(charsets_premasA(ps[i],as,ord))]),setify1([0]))
1114 )$
1115 charsets_premasA ( f,as,ord):=block( [ remd,i],
1116    remd : f,
1117    for i : length ( as)  thru 1 step  -  1 do
1118       remd : charsets_premA ( remd, as[i], charsets_lvar ( as[i], ord) ),
1119    if remd # 0 then (
1120       num( remd / lcoeff ( remd,ord) ))
1121    else ( 0)
1122 )$
1123 charsets_checkwith ( ps1,ps2):=block( [ rs,i,j,r,g,last],
1124    rs : charsets_setdifference ( ps1, ps2),
1125    if rs =  setify1([]) then ( true)
1126    else (
1127      if ps2 =  setify1([])  then ( false)
1128      else (
1129          rs : setify1(map(charsets_pfactor,map(sqfr,rs))),
1130          last:true,
1131          for i : 1 thru length ( rs)  do (
1132                r : rs[i],
1133                for j : 1 thru length ( ps2)  do (
1134                   g:gcd ( r, ps2[j]), r:r/g,
1135                   if ratnump(r) then ( return(last:true))),
1136                if  not  ratnump(r) then ( return ( last:false))),
1137          last))
1138 )$
1139 /* replace the power of factors of polys in as by 1 if any */
1140
1141 charsets_nopower ( as):=block(
1142    [ i],
1143       if  not   charsets_operatorp(as, charsets_ListOrSet)
1144           then (
1145             if  charsets_operatorp(as,"^")
1146                 then (
1147                   part( as, 1))
1148                else if  charsets_operatorp(as,"*") then (
1149                   apply("*",maplist(charsets_nopower, as)))
1150                else (
1151                   as))
1152
1153          else (
1154              maplist(charsets_nopower, as))
1155 )$
1156 /* macro to define the nearly char set -- a medial set */
1157 makeCharsetnFunction(theName,theSetType)::=buildq([theName,theSetType],
1158    theName( _ps,ord, [theRemsetFunction]):=block([ ps,cs,rs],
1159       if theRemsetFunction = [] then
1160           theRemsetFunction:charsets_remseta
1161       else theRemsetFunction:first(theRemsetFunction),
1162       if length ( _ps)  < 2 then ( _ps)
1163       else (
1164          ps:charsets_unigcd(_ps,ord),
1165          cs : theSetType ( ps, ord),
1166          if 0 < charsets_class ( cs[1], ord) then (
1167             rs : charsets_remseta ( charsets_minus ( ps, cs), cs, ord))
1168          else ( return ( cs)),
1169          if rs =  setify1([])  then ( cs)
1170          else ( theName(  append(rs, cs), ord)))
1171   ))$
1172
1173 /* the nearly char set -- a medial set */
1174 makeCharsetnFunction(charsets_charsetn,charsets_basset)$
1175 /* the nearly weak-char set -- a weak-medial set */
1176 makeCharsetnFunction(charsets_wcharsetn,charsets_wbasset)$
1177
1178 /* the nearly quasi-char set -- a quasi-medial set */
1179 makeCharsetnFunction(charsets_qcharsetn,charsets_qbasset)$
1180
1181 makeFcharsetnFunction(theName,theSetType)::=buildq([theName,theSetType],
1182    theName( _ps,ord,fset1,[theRemsetFunction]):=block(
1183       [ ps,cs,rs,ts,fset2,fset3,last,fset],
1184       if theRemsetFunction = [] then
1185           theRemsetFunction:charsets_remseta
1186       else theRemsetFunction:first(theRemsetFunction),
1187       if length ( _ps)  < 2 then ([_ps,fset1])
1188          else (
1189             ps : charsets_unigcd(_ps,ord),
1190             cs : theSetType( ps, ord),
1191             fset2 : [ fset1[1], charsets_union ( fset1[2], charsets_initialset1 ( cs, ord) )],
1192             if 0 < charsets_class ( cs[1], ord)  then (
1193                   rs:charsets_remseta ( charsets_minus ( ps, cs), cs, ord),
1194                   if length(rs) # 0 then (
1195                       ml2(rs,ts,charsets_removecont(rs,ord)),
1196                       ml2(rs,fset3, charsets_removefactor ( rs, ord, fset2)),
1197                       fset3 : [charsets_union(fset3[1],ts),fset3[2]])
1198                   else (
1199                       rs:[],
1200                       fset3:fset2))
1201                else (
1202                   return ( [cs,fset2]))),
1203             if rs =  [] then (
1204                   [cs,fset3])
1205                else (
1206                   theName( append(rs,cs), ord, fset3))
1207 ))$
1208 /* the modified nearly char set -- a modified medial set */
1209 makeFcharsetnFunction(charsets_fcharsetn,charsets_basset)$
1210
1211 /* the modified nearly weak-char set -- a modified weak-medial set */
1212 makeFcharsetnFunction(charsets_fwcharsetn,charsets_wbasset)$
1213
1214 /* the modified nearly quasi-char set -- a modified quasi-medial set */
1215 makeFcharsetnFunction(charsets_fqcharsetn,charsets_qbasset)$
1216
1217 /* the triangular set of poly set ps -- a quasi-medial set */
1218
1219 charsets_triset ( _ps,ord):=block(
1220    [ ps,i, flag,r],
1221       if length ( _ps)  < 2 then ( _ps)
1222          else (
1223             init_qs(ord),
1224             ps:charsets_unigcd(_ps,ord),
1225             for i : 0 thru length ( ord)  do charsets_qs[i] :  [],
1226             for i in ps do (
1227                charsets_qs[charsets_class ( i, ord) ] :
1228                   endcons(i,charsets_qs[charsets_class ( i, ord) ])),
1229             flag:false,
1230             for i : length ( ord)  thru 1 step  -  1 do
1231                if charsets_qs[0] =  []
1232                    then (
1233                      charsets_subtriset ( i, ord))
1234                   else ( return(flag:true)),
1235             if flag = true then return (  [ 1 ] ),
1236             if charsets_qs[0] #  []
1237                 then (
1238                    [ 1 ])
1239                else (
1240                    r:[],
1241                    for i:1 thru length(ord) do r:append(r,charsets_qs[i]),
1242                    r))
1243
1244 )$
1245 /* subroutine for triset */
1246
1247 charsets_subtriset ( i,ord):=block(
1248    [ ss,ss1,j,p],
1249       if 1 < length ( charsets_qs[i])
1250           then (
1251             ml2(ss1,ss, charsets_sort( charsets_qs[i], charsets_trank, ord)),
1252             charsets_qs[i] :  [ ss1 ],
1253             for j in ss do (
1254             p : charsets_prem ( j, ss1, ord[i]),
1255             if p # 0 then (
1256               charsets_qs[charsets_class ( p, ord) ] :
1257               endcons(num( p/ lcoeff ( p,ord)),
1258                    charsets_qs[charsets_class( p, ord)]))),
1259             charsets_subtriset ( i, ord))
1260 )$
1261 /* remove factor g from f until f has no factor g */
1262 /*      if g is removed then assign true to ja       */
1263
1264 charsets_movefactor ( f,g,ord):=block([ fg,ja,last_result],
1265    if not charsets_operatorp(f,charsets_ListOrSet) and not ratnump(g)
1266    and ((fg:first(last_result:divide(f,g,
1267            charsets_lvar(g,ord))),second(last_result)) = 0) then (
1268        ja:if 0 < charsets_class ( g, ord) then true else false,
1269        [first(charsets_movefactor ( fg, g, ord)),ja])
1270    else [f,false]
1271 )$
1272 /* remove possible factors in fset1[1] and fset1[2] from polys in ps */
1273 /*      where fset1[1] contains all factors removed before  */
1274 /*      if any poly in fset1[2] is removed, it is added to fset1[1] */
1275 /*      fset1 is assigned to fset at the end of the procedure */
1276
1277 charsets_removefactor ( ps,ord,fset1):=block(
1278    [ k,rr,ja,fset2,fs,qs,rs,r,fset],
1279    if  not   charsets_operatorp(ps,charsets_ListOrSet) then (
1280       qs :  setify1([ ps ]))
1281    else ( qs : ps),
1282    rs :  setify1([]),
1283    fs : fset1,
1284    for r in qs do (
1285       rr : r,
1286       fset2 :  setify1([]),
1287       for k in fs[1] do ml2(rr,ja,charsets_movefactor ( rr, k, ord)),
1288       for k in fs[2] do
1289          if rr # k then (
1290             ml2(rr,ja,charsets_movefactor ( rr, k, ord)),
1291             if ja then (
1292                fset2 : charsets_adjoin2(num( k / lcoeff( expand( k),ord)), fset2))),
1293       fs :  [ charsets_union ( fs[1], fset2), charsets_setdifference ( fs[2], fset2)  ],
1294       rs :  charsets_adjoin2( rr, rs)),
1295    rs : charsets_nopower ( rs),
1296    fset : fs,
1297    if  not charsets_operatorp(ps,charsets_ListOrSet) then ( [rs[1],fset])
1298    else ( [rs,fset])
1299 )$
1300 /* remove contents of all polys in ps wrt leading variables  */
1301 /*      the set of removed factors is assigned to ms */
1302
1303 charsets_removecont ( ps,ord):=block(
1304    [ qs,fs,i,cc,pp, ms],
1305      if  charsets_operatorp(ps, charsets_ListOrSet) then (
1306          if charsets_class ( ps[1], ord)  = 0 then [ps,setify1([])]
1307          else (
1308             qs :  [],
1309             fs :  setify1([]),
1310             for i : 1 thru length ( ps)  do (
1311                if charsets_class ( ps[i], ord)  > 0 then (
1312                ml2(cc,pp, content ( ps[i], charsets_lvar ( ps[i], ord))),
1313                cc:charsets_nopower(cc),
1314                if 0 < charsets_class ( cc, ord) then (fs : charsets_adjoin2(cc, fs)),
1315                pp :  expand( pp),
1316                qs :  endcons( pp,qs))),
1317             [qs,fs]))
1318      else (
1319         ml2(qs,ms,charsets_removecont([ps],ord)),
1320         [first(qs),ms])
1321 )$
1322 /* the modified triangular set of poly set ps -- a modified quasi-medial set */
1323
1324 charsets_ftriset ( ps,ord,fset1):=block([i,fset2,var,last],
1325    charsets_fact : 1,
1326    if length ( ps)  < 2 then ([ps, fset1])
1327    else (
1328       fset2 :  setify1([]),
1329       init_qs(ord),
1330       for i : 0 thru length ( ord)  do
1331                charsets_qs[i] :  [],
1332       last:[],
1333       for i in ps do (
1334          charsets_qs[charsets_class ( i, ord) ] :
1335              endcons( i, charsets_qs[charsets_class ( i, ord) ])),
1336       var : listofvars ( ps),
1337       for i : length ( ord)  thru 1 step  -  1 do
1338          if charsets_qs[0] =  [] then (
1339             last : charsets_fsubtriset ( i, ord, ps, var),
1340             fset2 : charsets_union ( fset2, last)
1341          )
1342          else ( return ( last:[[ 1 ],[ fset2,  setify1([]) ]])),
1343       if last # [] then return(last),
1344       if charsets_qs[0] #  [] then ([[1], [ fset2,  setify1([]) ]])
1345       else ([flatten_qs(1,length( ord)),[fset2,setify1([])]]))
1346 )$
1347 /* subroutine for ftriset */
1348
1349 charsets_fsubtriset ( i,ord,ps,var):=block(
1350    [ fset2,ss,ss1,j,k,l,p,pp,qq,ja,ts,last],
1351    if 1 < length ( charsets_qs[i]) then (
1352       ml2(ss1,ss,charsets_sort( charsets_qs[i], charsets_trank, ord)),
1353       charsets_qs[i] :  [ ss1 ],
1354       fset2 :  setify1([]),
1355       qq: setify1( maplist(lambda([p],num(p/lcoeff(expand(p),ord))),ps)),
1356       for j in ss do (
1357          pp : charsets_prem ( j, ss1, ord[i]),
1358          if pp # 0 then (
1359             if pp # charsets_fact and charsets_fsubtrisetsub(var,charsets_fact)
1360                 and  not  member(charsets_fact, qq) then (
1361                ml2(p,ts, charsets_removecont ( pp, ord)),
1362                fset2 :  setify1(append(fset2,ts)), /* could be charsets_union */
1363                ml2(p,ja,charsets_movefactor ( p, charsets_fact, ord)),
1364                if ja then (
1365                   fset2:charsets_adjoin2(num(charsets_fact
1366                       /lcoeff(expand(charsets_fact),ord)),fset2)
1367                ))
1368             else ( p : pp),
1369             for k in var do (
1370                if p # k and  not  member(k, qq) then (
1371                   ml2(p,ja,charsets_movefactor ( p, k, ord)),
1372                   if ja then (
1373                      fset2 :  charsets_adjoin2( k, fset2)))),
1374             charsets_qs[charsets_class ( p, ord) ] :
1375                endcons( charsets_nopower(num(p/lcoeff(p,ord))),
1376                charsets_qs[charsets_class ( p, ord)])
1377          )),
1378          last:charsets_initial ( ss1, ord),
1379          charsets_fact : num( last / lcoeff ( expand ( last),ord ) ),
1380          last : charsets_fsubtriset ( i, ord, ps, var),
1381          charsets_union ( fset2, last ))
1382       else ( setify1([]))
1383 )$
1384 /* subroutine for fsubtriset */
1385
1386 charsets_fsubtrisetsub ( aa,bb):=block(
1387    [ i,flag:true],
1388    for i : 1 thru length ( aa)  do
1389       if subst( aa[i] = 0, bb) = 0 then ( return ( flag:false)),
1390    flag
1391 )$
1392 /* reduce a triangular set into an ascending set */
1393
1394 charsets_trisetc ( ps,ord):=block(
1395    [ _ind,cs],
1396       charsets_cs : charsets_triset ( ps, ord),
1397       ml2(cs,_ind, charsets_subtrisetc ( ord,  setify1([]))),
1398       if _ind = 0 then (
1399          charsets_charsetn (  append(ps,cs), ord))
1400       else cs
1401 )$
1402 /* reduce a triangular set into an ascending set with factors moved */
1403
1404 charsets_ftrisetc ( ps,ord,fset1):=block( [ i,_ind,cs,fs,last,fset],
1405       ml2(charsets_cs,fs, charsets_ftriset ( ps, ord, fset1)),
1406       fset : fs,
1407       if fs[1] #  setify1([]) then (
1408              last:setify1(maplist(charsets_pfactor, fs[1])))
1409          else (
1410              last:setify1([])),
1411       ml2(cs,_ind, charsets_subtrisetc ( ord, last)),
1412       if _ind = 0 then ( charsets_fcharsetn ( append(ps,cs), ord, fset1))
1413       else ( [cs,fset])
1414 )$
1415 /* subroutine for ftrisetc */
1416
1417 charsets_subtrisetc ( ord,var):=block([ r,i,j,cs,_ind,nops],
1418    nops:if setify1(charsets_cs) = [[]] then 0 else length(charsets_cs),
1419    /* print("subtrisetc nops",nops), */
1420    if nops = 0 then ( cs :  [])
1421    else (cs :  [ charsets_cs[1] ]),
1422    if 1 < nops  then (
1423       for i : 2 thru nops  do (
1424          r : charsets_premas ( charsets_cs[i], cs, ord),
1425          if charsets_class( r, ord) # charsets_class( charsets_cs[i],ord) then (
1426             _ind : 0,
1427             if r # 0 then (
1428                cs : endcons( charsets_nopower( num(r/lcoeff(r,ord))),cs)),
1429                return(false))
1430             else (
1431                for j in var do r : first(charsets_movefactor ( r, j, ord)),
1432                cs : endcons( charsets_nopower( num( r / lcoeff ( r,ord))),cs),
1433                if charsets_class( r, ord) #
1434                       charsets_class( charsets_cs[i], ord) then (
1435                   return(false))))),
1436    [cs,_ind]
1437 )$
1438 /* the set of nonconstant initials of as   */
1439 /*    with certain repeated factors cancelled */
1440
1441 charsets_initialset1 ( as,ord):=block(
1442    [ i,__is,iss,last],
1443       __is :  setify1([]),
1444       for i in as do (
1445          last:charsets_initial ( i, ord),
1446          if 0 < (charsets_class ( last, ord)) then (
1447                __is : charsets_adjoin2(charsets_pfactor ( last), __is))),
1448       __is : charsets_compress ( __is, ord),
1449       iss :  setify1([]),
1450       for i in __is do
1451          if 0 < charsets_class ( i, ord)
1452              then (
1453                iss : charsets_adjoin2(i, iss)),
1454       iss
1455 )$
1456 /* compress some repeated factors */
1457
1458 charsets_compress ( ps,ord):=block([_is,is1,i,j,ss],
1459    _is : ps,
1460    if 1 < length ( _is) then (
1461       is1 :  [],
1462       for i : 1 thru length ( _is)  - 1 do (
1463          ss : _is[i],
1464          for j : i + 1 thru length ( _is)  do
1465             ss : first(charsets_movefactor ( ss, _is[j], ord)),
1466          if 0 < charsets_class ( ss, ord) then (
1467             is1 :  cons( charsets_pfactor( ss), is1))),
1468          is1 :  cons( _is[length ( _is) ], is1),
1469          _is :  setify1([]),
1470          if 1 < length ( is1) then (
1471             for i : 1 thru length ( is1)  - 1 do (
1472                ss : is1[i],
1473                for j : i + 1 thru length ( is1) do
1474                   ss : first(charsets_movefactor ( ss, is1[j], ord)),
1475                if 0 < charsets_class ( ss, ord) then (
1476                   _is :  setify1(charsets_adjoin2( charsets_pfactor ( ss),_is)))),
1477             setify1(charsets_adjoin2( is1[length ( is1) ],_is)))
1478          else ( setify1(is1)))
1479       else ( setify1(_is))
1480 )$
1481 /* the sequence of distinct factors of f */
1482
1483 charsets_pfactor ( f):=block(
1484    [ i,last],
1485       if integerp(f) then ( setify1([]))
1486          else if  charsets_operatorp(f,"^") then (
1487             last:part( f, 1),
1488             num( last / lcoeff ( last) ))
1489          else if  charsets_operatorp(f,"*") then (
1490             maplist(charsets_pfactor,f))
1491          else (
1492             num( f / lcoeff ( expand ( f) ) ))
1493 )$
1494 /* the sequence of factors of f */
1495
1496 charsets_sfactor ( f):=block(
1497    [ i,last],
1498       if integerp(f) then (setify1([]))
1499          else if  charsets_operatorp(f,"^") then (
1500             last:part( f, 1),
1501             makelist(numer(last/lcoeff(last)),i,1,part( f, 2)))
1502          else if  charsets_operatorp(f,"*") then (
1503             map(charsets_sfactor,f))
1504          else (
1505             num( f / lcoeff ( f) ))
1506 )$
1507 /* the set of all nonconstant factors of initials of polys in as */
1508
1509 charsets_initialset ( as,ord):=block(
1510    [ i,__is,iss,last],
1511    __is :  setify1([]),
1512    for i in as do (
1513       last:charsets_initial ( i, ord),
1514       if 0 < (charsets_class ( last, ord)) then (
1515          __is :  setify1(cons( last, __is)))),
1516    __is : charsets_factorps ( __is),
1517    iss :  setify1([]),
1518    for i in __is do (
1519       if 0 < charsets_class ( i, ord) then (
1520          iss :  setify1(cons( i, iss)))),
1521    iss
1522 )$
1523 /* all irreducible nonconstant factors of a set of polynomials */
1524
1525 charsets_factorps ( ps):=block([ qs,i,j,q,c],
1526    qs :  setify1([]),
1527    for i in ps do (
1528       q : factor ( i),
1529       c : num(numfactor(q)),
1530       q : q/c,
1531       if (charsets_operatorp(q,"*") = true) or (charsets_operatorp(q,"/") = true) then (
1532          for j : 1 thru length ( q)  do (
1533             if  not  integerp(part( q, j)) then (
1534                if  charsets_operatorp(part( q, j),"^") then (
1535                   qs :  setify1(cons( c*num( part( part( q, j), 1)
1536                                   / lcoeff ( part( part( q, j), 1)) ),qs)))
1537                else (
1538                   qs :  setify1(cons( c*num( part( q, j)
1539                                   / lcoeff ( part( q, j)) ),qs))))))
1540       else if  charsets_operatorp(q,"^") then (
1541          qs :  setify1(cons(c*num( q/ lcoeff ( q)),qs)))
1542          else (
1543             if not integerp(q) then (
1544                qs : setify1(cons(c*num( q / lcoeff ( q) ),qs))))),
1545       qs
1546 )$
1547 /*#### Part II. Routines for Various Decompositions and Others ##### */
1548 /* the char series of poly set ps */
1549
1550 charsets_hack1(iss):=block([l1:first(iss),l2:second(iss)],
1551    [[l1[2],l1[3],l1[4],l1[1]],
1552     [l2[2],l2[3],l2[1],l2[4]]])$
1553
1554 charsets_hack2(a):=[second(a),first(a)]$
1555
1556 charsets_charseries ( ps,ord,medset):=block(
1557    [ qs,cs,iss,n,qsi,qhi,csno,ppi,qqi],
1558    if charsets_operatorp(ps[1], charsets_ListOrSet) then ( qhi : setify1(ps))
1559    else ( qhi :  setify1([setify1( ps)])),
1560    qsi :  setify1([]),
1561    csno : 0,
1562    ppi :  setify1([]),
1563    qqi :  setify1([]),
1564    for n : 0 while qhi #  setify1([]) do (
1565       qhi : sort (qhi, charsets_lenord),
1566       qs : qhi[1],
1567       ppi : charsets_select ( ppi, length ( qs) ),
1568       qqi :  setify1(append(ppi[2], qqi)),
1569       if n = 0 then ( ppi :  setify1([]))
1570       else ( ppi :  charsets_adjoin2(qs,ppi[1])),
1571       cs : charsets_charseta ( qs, ord, medset),
1572       if 1 < charsets_printlevel then (
1573          csno : csno + 1,
1574          print("Characteristic set produced in charseries",
1575              csno,length ( qhi),length ( qsi),length ( qs)),
1576          print ( cs)),
1577       if 0 < charsets_class ( cs[1], ord) then (
1578          iss : charsets_initialset ( cs, ord),
1579          if charsets_simpa ( iss, cs, ord)  # 0 then (
1580             qsi : charsets_adjoin2( cs, qsi)),
1581          iss : charsets_adjoin3 ( iss, qs, qqi))
1582       else ( iss :  setify1([])),
1583       /* if read("reorder") = 1 then iss:charsets_hack1iss, */
1584       if 1 < length ( qhi) then (
1585          qhi :  setify1(append(iss, rest(qhi))))
1586       else ( qhi : iss)),
1587    if qsi =  setify1([]) then ( [])
1588    else ( sort ( charsets_contract ( qsi, ord, 0), charsets_lenord))
1589 )$
1590 /* the char series of poly set ps -- allowing to remove factors */
1591
1592 charsets_fcharser ( ps,ord,medset):=block(
1593    [ qs,cs,iss,n,qhi,qsi,factorset,csno,ppi,qqi,last],
1594    if  charsets_operatorp(ps[1], charsets_ListOrSet) then ( qhi :  setify1(ps))
1595    else ( qhi :  setify1([setify1(ps)])),
1596    qsi :  setify1([]),
1597    csno : 0,
1598    ppi :  setify1([]),
1599    qqi :  setify1([]),
1600    for n : 0 while qhi #  setify1([]) do (
1601       qhi : sort (qhi, charsets_lenord),
1602       qs : qhi[1],
1603       ppi : charsets_select ( ppi, length ( qs) ),
1604       qqi :  setify1(append(ppi[2], qqi)),
1605       if n = 0 then ( ppi :  setify1([]))
1606       else ( ppi :  charsets_adjoin2(qs,ppi[1])),
1607       if length ( qs)  - 3 < length ( ord) then (
1608             cs : charsets_fcharseta (qs, ord, medset),
1609             factorset : cs[2],
1610             cs : cs[1])
1611       else (
1612             cs:charsets_charseta (qs, ord, medset),
1613             factorset:setify1([])),
1614       if 1 < charsets_printlevel then (
1615             csno : csno + 1,
1616             print("Characteristic set produced in fcharser",
1617                csno,length ( qhi),length ( qsi),length ( qs)),
1618             print ( cs)),
1619       if 0 < charsets_class ( cs[1], ord)  then (
1620             iss : charsets_initialset ( cs, ord),
1621             if charsets_simpa ( iss, cs, ord)  # 0 then (
1622                qsi : charsets_adjoin2( cs, qsi)),
1623             iss : charsets_union ( iss, charsets_factorps ( factorset) ))
1624          else (
1625             iss : charsets_factorps ( factorset)),
1626       iss : charsets_adjoin3 ( iss, qs, qqi),
1627       if 1 < length ( qhi)  then (
1628          qhi :  setify1(append(iss, rest(qhi))))
1629          else ( qhi : iss)),
1630    if qsi =  setify1([])  then ( [])
1631    else ( sort ( charsets_contract ( qsi, ord, 0), charsets_lenord))
1632 )$
1633 /* the extended char series of poly set ps */
1634
1635 charsets_excharser ( ps,ord,medset):=block(
1636    [qs,cs,_is,iss,i,j,qsi,qhi,r,rr,last],
1637    if  charsets_operatorp(ps[1], charsets_ListOrSet) then ( qhi :  setify1([ ps ]))
1638    else ( qhi :  setify1([  [ ps, 1 ]  ])),
1639    qsi :  setify1([]),
1640    while (qhi #  setify1([])) do (
1641       qs : qhi[1][1],
1642       cs : charsets_charseta (qs, ord, medset),
1643       if 1 < charsets_printlevel then (
1644                print("Characteristic set produced in excharser"),
1645                print ( cs)),
1646       if 0 < charsets_class ( cs[1], ord) then (
1647          _is : charsets_initialset ( cs, ord),
1648          rr : charsets_nopower( charsets_prod( charsets_adjoin2(qhi[1][2],_is))),
1649          last:charsets_premas ( rr, cs, ord),
1650          last:factor(last),
1651          r : charsets_simp ( last, cs, ord),
1652          if r # 0 then (
1653             if r = 1 then ( qsi : setify1(cons( cs, qsi )))
1654             else ( qsi : charsets_adjoin2([ cs, charsets_simpb( r, rr)], qsi))))
1655       else ( _is : []),
1656       iss :  setify1([]),
1657       if length ( ord)  <= length ( ps)  + 1 then (
1658          for i in _is do ( iss : setify1(cons([ charsets_adjoin2( i,qs),qhi[1][2]],iss))))
1659       else (
1660          for i : 1 thru length ( _is)  do (
1661             last:if i = 1 then 1 else product ( _is[j], j, 1, i - 1),
1662             iss:charsets_adjoin2([charsets_adjoin2(_is[i],qs),last*qhi[1][2]],iss)
1663          )
1664       ),
1665       if 1 < length ( qhi) then (
1666          qhi : charsets_union(iss, rest(qhi)))
1667       else ( qhi : iss)
1668    ),
1669    if qsi =  setify1([]) then ( [])
1670    else ( sort ( qsi, charsets_lenord))
1671 )$
1672 /* simplify r to rr so that Zero(cs/r) = Zero(cs/rr) holds still  */
1673 charsets_simp ( r,cs,ord):=block([rr,i,fs,j,_ind,last],
1674    if r = 0 then ( 0)
1675    else (
1676       fs : setify1(charsets_pfactor( r)),
1677       rr : 1,
1678       for i in fs do (
1679          if 0 < charsets_class ( i, ord)  then (
1680                _ind : 0,
1681                for j in cs do (
1682                   last:subst(i = 0,j),
1683                   if last = 0 then (
1684                      if charsets_class(first(charsets_movefactor(j,i,ord)),
1685                             ord)  = 0 then (
1686                         _ind :  -  1,
1687                         return(true)))
1688                      else if charsets_class ( last, ord)  = 0 then (
1689                         _ind : 1,
1690                         return(true))),
1691                if _ind = 0 then ( rr : rr * i)
1692                else if _ind =  -  1 then (
1693                   return(true)))),
1694       if _ind =  -  1 then ( 0) else ( rr))
1695 )$
1696 /* check whether Zero(cs/fs) is empty */
1697
1698 charsets_simpa ( fs,cs,ord):=block([ i,j,ds,flag],
1699    if length ( cs)  = 1 then ( 1)
1700    else (
1701       flag:1,
1702       ds : makelist(cs[i],i,1,length ( cs)  - 1),
1703       for i in fs do (
1704          for j in ds do (
1705             if subst(0, i, j) = 0 then (
1706                if charsets_class( first(charsets_movefactor(j,i,ord)),ord) = 0
1707                   then ( return ( flag:0)))),
1708          if(flag = 0)then return(flag)),
1709       flag)
1710 )$
1711 /* the simpler one of a and b */
1712
1713 charsets_simpb ( a,b):=block(
1714       if charsets_measure ( a)  < charsets_measure ( b) then ( a) else ( b)
1715 )$
1716 /* the measure of complexity of a according to number of terms */
1717
1718 charsets_measure ( a):=block( [ i],
1719    if  charsets_operatorp(a,"^") then ( length ( part( a, 1)))
1720    else if  charsets_operatorp(a,"*") then (
1721       sum ( charsets_measure ( part( a, i)), i, 1, length ( a) ))
1722    else if variablep(a) then 1
1723    else ( length ( a))
1724 )$
1725 /* the extended char series of poly set ps -- allowing to remove factors */
1726
1727 charsets_fexcharser ( ps,ord,medset):=block(
1728    [ qs,cs,_is,iss,n,i,j,qhi,qsi,r,rr,factorset,last],
1729    if  charsets_operatorp(ps[1], charsets_ListOrSet) then ( qhi :  setify1([ ps ]))
1730    else ( qhi :  setify1([  [ ps, 1 ]  ])),
1731    qsi :  setify1([]),
1732    for n : 0 while qhi #  setify1([])  do (
1733       qs : qhi[1][1],
1734       if n = 0 then ml2(cs,factorset,charsets_fcharseta(qs, ord, medset))
1735       else (
1736          cs:charsets_charseta (qs, ord, medset),
1737          factorset:setify1([])),
1738       if 1 < charsets_printlevel then (
1739          print("Characteristic set produced"),
1740          print ( cs)),
1741       if 0 < charsets_class ( cs[1], ord) then (
1742          _is : charsets_union(charsets_initialset(cs,ord),charsets_factorps(factorset)),
1743          rr : charsets_nopower( charsets_prod( charsets_adjoin2(qhi[1][2],_is))),
1744          last:charsets_premas ( rr, cs, ord),
1745          last:factor(last),
1746          r : charsets_simp ( last, cs, ord),
1747          if r # 0 then (
1748             if r = 1 then ( qsi :  setify1(cons( cs, qsi )))
1749             else ( qsi : setify1(endcons([ cs, charsets_simpb( r, rr)],qsi)))))
1750       else ( _is : charsets_factorps ( factorset)),
1751       iss :  setify1([]),
1752       if length ( ord)  <= length ( ps)  + 1 then (
1753          for i in _is do ( iss : setify1(endcons([charsets_adjoin2( i,qs),qhi[1][2]],iss))))
1754       else (
1755          for i : 1 thru length ( _is)  do (
1756             last:if i = 1 then 1 else product ( _is[j], j, 1, i - 1),
1757             iss:charsets_adjoin2([charsets_adjoin2(_is[i],qs),last*qhi[1][2]],iss)
1758          )
1759       ),
1760       if 1 < length ( qhi) then qhi : charsets_union(iss,rest(qhi))
1761       else qhi : iss
1762    ),
1763    if qsi =  setify1([]) then ( [])
1764    else (sort ( qsi, charsets_lenord))
1765 )$
1766 /* the irreducible char series of poly set ps */
1767 /*     using new factorization method if m=1, Hu-Wang's method if m=-1 */
1768 /*     and normalized char set if m=0     */
1769
1770 charsets_irrcharser ( ps,ord,medset,[m]):=block(
1771    [ qs,cs,cst,_is,iss,n,ts,qsi,qhi,pi,factorset,ppi,qqi,csno,_ind,mind,fset],
1772    if length(m) = 0 then ( _ind : 1) else ( _ind : first(m)),
1773    if  charsets_operatorp(ps[1], charsets_ListOrSet) then ( qhi :  setify1([ (ps) ]))
1774    else ( qhi :  setify1([ ps ])),
1775    qsi :  setify1([]),
1776    pi :  setify1([]),
1777    csno : 0,
1778    ppi :  setify1([]),
1779    qqi :  setify1([]),
1780    if medset = charsets_basset then ( mind : true) else ( mind : false),
1781    for n : 0 while qhi #  setify1([])  do (
1782       qhi : sort ( qhi, charsets_nopsord),
1783       qs : qhi[1],
1784       ppi : charsets_select ( ppi, length ( qs) ),
1785       qqi :  charsets_union(ppi[2], qqi),
1786       if n = 0 then ( ppi :  setify1([]))
1787       else ( ppi : charsets_adjoin2(qs, ppi[1])),
1788       if length ( qs)  - 3 < length ( ord) then (
1789          if not mind then (
1790             ml2(cs,fset, apply( assoc(medset,charsets_fcharsetsub_alist),
1791                [ qs, ord, [ setify1([]), listofvars( qs)]])),
1792             factorset:first(fset))
1793          else (
1794             ml2(cs,factorset, charsets_fcharseta ( qs, ord, medset)),
1795             if 1 < charsets_printlevel and _ind = 1 then (
1796                csno : csno + 1,
1797                print("Characteristic set produced in irrcharser",csno,
1798                   length ( qhi),length ( qsi),length ( qs)),
1799                print ( cs))),
1800          if _ind = 0 then (
1801             cs :  [ charsets_fcnormal ( cs, ord, 2)  ],
1802             if 1 < length ( cs) then (
1803                factorset : charsets_union ( factorset, part( cs[2], 2))),
1804             cs : cs[1]))
1805       else (
1806          if mind then ml2(cs,factorset,
1807             charsets_removecont ( charsets_charseta ( qs, ord, medset), ord))
1808          else ml2(cs,factorset,charsets_removecont( medset (qs,ord),ord))),
1809       if 0 < charsets_class ( cs[1], ord) then (
1810          ts:charsets_irras ( cs, ord, _ind),
1811          if ts[2] = 0 then (
1812             if  not  mind then (
1813                if  medset = charsets_autored then
1814                   cs:charsets_charseta(charsets_union(qs,cs),ord, charsets_charsetn)
1815                else (
1816                   if  not charsets_subset ( cs, qs) then
1817                      cs : charsets_charseta( charsets_union(qs, cs), ord, medset)),
1818                if 1 < charsets_printlevel and _ind = 1 then (
1819                   csno : csno + 1,
1820                   print("Characteristic set produced in irrcharser",
1821                       csno,length ( qhi),length ( qsi),length ( qs)),
1822                   print ( cs))),
1823             if  not  member(cs, pi) then (
1824                pi :  charsets_adjoin2(cs, pi),
1825                if 0 < charsets_class ( cs[1], ord) then (
1826                   ts:charsets_irras ( cs, ord, _ind),
1827                   if ts[2] = 0 then (
1828                      qsi : charsets_adjoin2( cs, qsi),
1829                      if length ( cs)  = length ( ord) then (
1830                         _is : charsets_factorps ( factorset))
1831                      else (
1832                         _is : charsets_union ( charsets_initialset ( cs, ord),
1833                            charsets_factorps ( factorset) )),
1834                      iss : charsets_adjoin3 ( _is, qs, qqi)))
1835                else (
1836                   iss : charsets_adjoin3(charsets_factorps(factorset),qs,qqi)))
1837             else (iss:charsets_adjoin3(charsets_factorps(factorset), qs, qqi))),
1838          if ts[2] # 0 then ( _is : charsets_factorps ( factorset),
1839             if 1 < ts[2] then (
1840                cst : part(cs, makelist(i,i,1,ts[2] - 1)),
1841                _is : charsets_union ( _is, charsets_initialset ( cst, ord) ),
1842                iss : charsets_union ( charsets_adjoin3 ( _is, qs, qqi),
1843                   charsets_adjoinb ( ts[1], qs, qqi, cst) ))
1844             else ( iss : charsets_adjoin3 ( charsets_union ( _is, ts[1]), qs, qqi))))
1845       else ( iss : charsets_adjoin3 ( charsets_factorps ( factorset), qs, qqi)),
1846       if 1 < length ( qhi) then ( qhi :  charsets_union(iss, rest(qhi)))
1847       else ( qhi : iss),
1848       last:if qsi =  setify1([]) then ([])
1849       else (sort( charsets_contract ( qsi, ord, 1), charsets_lenord))),
1850   last
1851 )$
1852 /* test whether ps is a charsets_subset or sublist of qs */
1853
1854 charsets_subset ( ps,qs):=block( [ p,flag:true],
1855    for p in ps do
1856       if not member(p, qs) then ( return ( flag:false)),
1857    flag
1858 )$
1859 /* subroutine for irrcharser, qirrcharser and others */
1860
1861 charsets_adjoinb ( _is,qs,qh,cs):=block( [ iss,i,j,_ind,qhi,itt,qs_],
1862    iss :  setify1([]),
1863    qhi : charsets_setdifference ( qh, (qs_:setify1(qs))),
1864    if _is #  setify1([]) then (
1865       for i in _is do (
1866          itt :  setify1(cons( i, charsets_union(qs_, setify1(cs)))),
1867          _ind : 0,
1868          if 0 < length ( qhi) then (
1869             for j in qhi while _ind = 0 do (
1870                if charsets_subset ( j, itt) then ( _ind : 1))),
1871          if _ind = 0 then ( iss : charsets_adjoin2(itt,iss)))),
1872    iss
1873 )$
1874 /* examine the irreducibility of as for irrcharser */
1875 /* This function should be memoized */
1876 charsets_irras ( as,ord,inda):=block([ _ind,i,j,p,qs,n,fs,ja,dd,last],
1877    _ind : 1,
1878    ja : 0,
1879    dd : 1,
1880    for i : 1 thru length ( as)  do (
1881       p : factor ( as[i]),
1882       fs:if mapatom(p) then [p]
1883       else (
1884          if "-" = op(p) then -charsets_dfactors(-p) else charsets_dfactors(p)),
1885    /* print("irras after dfactors",fs), */
1886       qs :  setify1([]),
1887       for j : 1 thru length ( fs)  do (
1888             if 0 < charsets_class ( fs[j], ord) then qs : charsets_adjoin2(fs[j],qs)),
1889    /* print("irras after qs loop",qs), */
1890       last:charsets_lvar ( p, ord),
1891       if charsets_degree ( qs[1], last)  < charsets_degree ( p, last) then (
1892    /* print("irras degree compare hit",qs[1],p), */
1893                ja : 1,
1894                _ind : 0,
1895                return(false))),
1896    if _ind = 1 and 1 < length ( as) then (
1897       for nn:(n:1)
1898          while (n:nn) < length(as) and charsets_degreel(as[n:nn],ord)=1 do (0),
1899       if n < length(as) then ml3(qs,ja,dd,charsets_irrassub( as, n, ord, inda))
1900       else ( ja : 0)),
1901    [ qs, ja,dd ]
1902 )$
1903 /* subroutine for irras */
1904
1905 charsets_irrassub ( as,n,ord,_ind):=block(
1906    [ m,qs,i,vv,dd,mm,ja],
1907    for mm : (m:n + 1) while (m:mm) <= length ( as)
1908          and charsets_degreel ( as[m:mm], ord)  = 1 do (0),
1909    dd : 1,
1910    if m <= length ( as) then (
1911       vv : charsets_lvar ( as[m], ord),
1912       if _ind = -1 then (
1913          qs : charsets_factoras ( as[m],  makelist(as[i],i,1,m - 1), ord))
1914       else (
1915          qs : charsets_cfactor ( as[m],  makelist(as[i],i,1,m - 1), ord),
1916          dd : denom ( qs),
1917          qs :  setify1(charsets_qfactor ( num( qs), ord))),
1918       if max ( map(lambda([q],charsets_degree (q, vv)),qs))
1919                 = charsets_degree ( as[m], vv) then (
1920          ml3(qs,ja,dd, charsets_irrassub ( as, m, ord, _ind)))
1921       else ( ja : m))
1922    else ( ja : 0),
1923    [qs,ja,dd]
1924 )$
1925 /* collect distinct nonconstant factors of a polynomial q */
1926
1927 charsets_dfactors ( q):=block( [ qs,j],
1928    if  charsets_operatorp(q,"*") then (
1929       qs :  setify1([]),
1930       for j : 1 thru length ( q)  do
1931          if  not  integerp(part( q, j)) then (
1932             if  charsets_operatorp(part( q, j),"^") then (
1933                qs :  setify1([ qs, num( part( part( q, j), 1)
1934                       / lcoeff ( part( part( q, j), 1)))]))
1935             else ( qs : charsets_adjoin2( num( part( q, j)/lcoeff( part( q, j))),qs))))
1936    else if  charsets_operatorp(q,"^") then (
1937       qs :  setify1([ num( part( q, 1) / lcoeff ( part( q, 1)) )  ]))
1938    else (
1939       if  not  integerp(q) then (
1940          qs : setify1([ num( q / lcoeff ( q))]))),
1941    qs
1942 )$
1943 /* normalize ascending set cs wrt ord */
1944
1945 charsets_fcnorm ( cs,ord,m):=block(
1946    [ as,i],
1947       as :  [ cs[1] ],
1948       for i : 2 thru length ( cs)  do
1949          as : charsets_fcnormal (  [ (as), cs[i] ], ord, m),
1950          if 1 < length (  [ as ] ) then ( return ( cs)) ,
1951       as
1952 )$
1953 /* normalize the last pol in an ascending set cs wrt ord */
1954
1955 charsets_fcnormal ( cs,ord,m):=block(
1956    [ n,ini,i,j,ggg,gg,ff,ccs,dd,cd,fs,nt,last],
1957    n : length ( cs),
1958    if n < 2 then ( cs)
1959    else (
1960       dd : cs[n],
1961       nt : length ( expand ( dd) ),
1962       for i : n - 1 thru 1 step  -  1 do (
1963          ini : charsets_initial ( dd, ord),
1964          if 0 < charsets_degree ( ini, charsets_lvar ( cs[i], ord) ) then (
1965             ggg : charsets_gcdex ( cs[i], ini, charsets_lvar ( cs[i], ord) ),
1966             gg : ggg[3],
1967             if 0 < charsets_degree ( gg, charsets_lvar ( cs[i], ord) ) then (
1968                ff : cs[i],
1969                gg :  setify1([ charsets_pfactor ( gg)  ]),
1970                cd :  setify1([]),
1971                for j : 1 thru length ( gg)  do
1972                   if charsets_class( gg[j], ord) = charsets_class ( cs[i], ord)
1973                      then (
1974                      ff : charsets_nopower( first(
1975                         charsets_movefactor( ff, gg[j], ord)) ),
1976                      cd :  setify1([ (cd), gg[j] ])) ,
1977                if charsets_class ( ff, ord)  = 0 then (
1978                   ccs :  [  [ 1 ]  ])
1979                else (
1980                   cs[i]:ff,
1981                   ccs : [charsets_fcnormal( cs, ord, m)]),
1982                if length ( ccs)  = 1 then (
1983                   return ([ccs[1], "common divisors", cd]))
1984                else (
1985                   return ([ccs[1], "common divisors",
1986                      setify1([cd, (part( ccs[2], 2)) ]) ])))
1987             else (
1988                dd : charsets_prem( dd*ggg[2], cs[i], charsets_lvar(cs[i],ord)),
1989                dd :first(charsets_movefactor(dd,
1990                   charsets_initial(cs[i],ord),ord)),
1991                dd : charsets_nopower ( dd),
1992                if 4 * m * nt < length ( expand ( dd) ) then ( return ( cs))))),
1993       ccs : makelist(cs[i],i,1,n - 1),
1994       block([temp], ml2(temp,dd,content(dd,charsets_lvar(dd,ord))),last:temp),
1995       fs : setify1([charsets_pfactor(last)]),
1996       gg :  setify1([]),
1997       for i : 1 thru length ( fs)  do
1998          if 0 < charsets_class (fs[i],ord) then ( gg : setify1([gg, fs[i]])),
1999       gg : charsets_prod ( gg),
2000       ini : charsets_initialset ( ccs, ord),
2001       for i : 1 thru length ( ini)  do
2002          gg : first(charsets_movefactor ( gg, ini[i], ord)),
2003       dd : charsets_nopower ( gg)  * dd,
2004       gg :  [ charsets_pfactor ( num( dd / lcoeff ( expand ( dd),ord ) ) )  ],
2005       dd : 1,
2006       for i : 1 thru length ( gg)  do
2007          if 0 < charsets_class ( gg[i], ord) then ( dd : dd * gg[i]) ,
2008       if m * nt < length ( expand ( dd) ) then (
2009          if 1 < charsets_printlevel then (
2010             print("normalization fails:",nt,length ( expand ( dd) ))) ,
2011          cs)
2012       else ( endcons(dd,ccs)))
2013 )$
2014 /* the modified gcdex for fcnormal */
2015 /* This function was memoized, however maybe it can be replaced by the maxima gcd */
2016 charsets_gcdex ( A,B,x):=block(
2017 [ m,pm,cc,cd,c,c1,c2,d,d1,d2,r,r1,r2,q,II,temp],
2018    if A = 0 then ( return (  [ 0, 1, B ] )) ,
2019    if B = 0 then ( return (  [ 1, 0, A ] )) ,
2020       ml2(cc,c,content( A, x)),
2021       ml2(cd,d,content( B, x)),
2022       /* this is a guess about how gcd/degree should work */
2023       II:block([g],g:gcd( c:expand ( c), d:expand ( d), x), c:c/g, d:d/g,
2024          if 1=g then setify1([]) else setify1([g])),
2025       pm : 1,
2026       c1 : 1,
2027       c2 : 0,
2028       d1 : 0,
2029       d2 : 1,
2030       while (d # 0) do (
2031          ml3(q,r,m,prem ( c, d, x)),
2032          r:first(divide ( r, pm)),
2033          r1:first(divide ( m * c1 - q * d1, pm)),
2034          r2:first(divide ( c2 * m - q * d2, pm)),
2035          c : d,
2036          c1 : d1,
2037          c2 : d2,
2038          d : r,
2039          d1 : r1,
2040          d2 : r2,
2041          pm : m
2042      ),
2043      map(expand,subst(II, [ c1 * cd, c2 * cc, c * cc * cd ]))
2044 )$
2045 /* the extended irreducible char series of poly set ps */
2046
2047 charsets_exirrcharser ( ps,ord,medset):=block(
2048    [qs,cs,_is,iss,n,i,j,qhi,qsi,r,rr,factorset,mind,fset,_ind,ts,den,last,temp],
2049    qhi:if charsets_operatorp(ps[1], charsets_ListOrSet) then setify1([ps])
2050    else setify1([[ps,1]]),
2051    mind:if medset = charsets_basset then true else false,
2052    qsi :  setify1([]),
2053    for n : 0 while qhi #  setify1([]) do (
2054       qs : qhi[1][1],
2055       if not mind then (
2056          if n < 20 then (
2057             ml2(cs,fset, apply( assoc(medset,charsets_fcharsetsub_alist),
2058                [ qs, ord, [ setify1([]), listofvars( qs)]])),
2059             factorset : fset[1])
2060          else (
2061             last:first(apply( assoc(medset,charsets_fcharsetsub_alist),
2062                [ qs, ord,[setify1([]),listofvars(qs)]])),
2063             ml2(cs,factorset, charsets_removecont ( last, ord))))
2064       else (
2065          if n < 20 then (
2066             ml2(cs,factorset, charsets_fcharseta( qs, ord, medset)))
2067          else (
2068             temp:charsets_charseta( qs, ord, medset),
2069             ml2(cs,factorset, charsets_removecont ( temp, ord))),
2070          if 1 < charsets_printlevel then (
2071             print("Characteristic set produced"),
2072             print ( cs))),
2073       if 0 < charsets_class ( cs[1], ord) then (
2074          /* in the original code ts has 2 elements, here it has 3, so
2075             confirm that the 3rd element is not used, if it is, then
2076             delete it after the assignment to den */
2077          ts:charsets_irras ( cs, ord, _ind),
2078          den:ts[3],
2079          if ts[2] = 0 then (
2080             if not mind then (
2081                if medset = charsets_autored then (
2082                   cs : charsets_charseta( setify1(append(qs,cs), ord, charsets_charsetn)))
2083                else (
2084                   if not  charsets_subset( cs, qs) then (
2085                      cs:charsets_charseta(setify1(append(qs,cs)),ord,medset))),
2086                 if 1 < charsets_printlevel then (
2087                         print("Characteristic set produced"),
2088                         print ( cs))),
2089             if 0 < charsets_class ( cs[1], ord) then (
2090                ts:charsets_irras ( cs, ord, _ind),
2091                den:ts[3],
2092                if ts[2] = 0 then (
2093                   _is : charsets_union( charsets_initialset( cs, ord), charsets_factorps( factorset)),
2094                   if length ( cs)  = length ( ord)  then (
2095                      rr : charsets_nopower ( qhi[1][2]))
2096                   else (
2097                      rr:charsets_nopower(charsets_prod(charsets_adjoin2(qhi[1][2],_is)))),
2098                      last:charsets_premas ( rr, cs, ord),
2099                      r : charsets_simp ( last, cs, ord),
2100                      if r # 0 then (
2101                            if r = 1 then (
2102                               qsi : charsets_adjoin2( cs, qsi))
2103                            else (
2104                               qsi : charsets_adjoin2([cs,charsets_simpb ( r, rr)],qsi))))
2105                         else (
2106                            _is : charsets_factorps ( factorset),
2107                            ts :  [ 1, 0 ]))),
2108          if ts[2] # 0 then (
2109             if 1 < ts[2] then (
2110                _is : charsets_union( charsets_initialset( setify1(
2111                    makelist(cs[i],i,1,ts[2]-1)),ord),
2112                    charsets_factorps(factorset)))
2113             else ( _is : charsets_factorps ( factorset))))
2114       else (
2115          _is : charsets_factorps ( factorset),
2116          ts :  [ 1, 0 ]),
2117       iss :  setify1([]),
2118       if length ( ord)  <= length ( ps)  + 1 then (
2119          for i in _is do (
2120             iss :  charsets_adjoin2( [charsets_adjoin2( i,qs ), qhi[1][2]],iss)))
2121       else (
2122          for i : 1 thru length ( _is)  do (
2123             last:if i = 1 then ( 1) else product ( _is[j], j, 1, i - 1),
2124             iss :  charsets_adjoin2([charsets_adjoin2(_is[i],qs), last * second(first(qhi))],iss))),
2125       if ts[2] # 0 and ts[1] #  setify1([])  then (
2126          if  not  mind then (
2127             last:"jello",
2128             last:if  medset =  charsets_autored then (
2129                charsets_charseta(setify1(append(qs,cs)),ord,charsets_charsetn))
2130             else (
2131                if  not  charsets_subset ( cs, qs)  then
2132                   charsets_charseta(setify1(append(qs,cs)), ord, medset)
2133                else cs),
2134             if last # cs then (
2135                if ts[2] = 1 then cs : qs
2136                else cs : setify1(append( qs, makelist(cs[i],i,1,ts[2] - 1))))
2137          ),
2138          for i in ts[1] do (
2139             iss : charsets_adjoin2([charsets_adjoin2(i,cs),
2140                charsets_prod(charsets_adjoin2(den,charsets_adjoin2(second(first(qhi)),_is)))],iss)),
2141             if 0 < charsets_class ( den, ord) and
2142                   ts[2] < charsets_class ( ts[1][1], ord)  then (
2143                iss :  charsets_adjoin2([ charsets_adjoin2( den, cs ), second(first(qhi))],iss))),
2144       if 1 < length ( qhi)  then ( qhi :  charsets_union(iss, rest(qhi)))
2145       else ( qhi : iss)),
2146    if qsi =  setify1([]) then ( [])
2147    else ( (sort (qsi,charsets_lenord)))
2148 )$
2149 /* subroutine for irrcharser, qirrcharser and others */
2150
2151 /* In maple the condition is n <= nops(i) and nops gives the same thing as maxima's */
2152 /* length for sums and products, excepts that for a single variable nops(x) = 1 while */
2153 /* length (x) produces an error. */
2154
2155 charsets_select ( ppi,n):=block([ i,pp,qq,lli],
2156    lli : if (i=listofvars(i)[1]) then 1 else length(i),
2157    pp :  setify1([]),
2158    qq :  setify1([]),
2159    for i in ppi do (
2160       if n <= lli then ( qq :  charsets_adjoin2( i, qq))
2161       else ( pp :  charsets_adjoin2( i, pp))),
2162    [ pp, qq ]
2163 )$
2164 /* subroutine for irrcharser, qirrcharser and others */
2165 /* name this function charsets_adjoin3 to distinguish from 2-arg
2166  * function in charsets_set.lisp.
2167  */
2168 charsets_adjoin3 ( _is,qs,qh):=block( [ iss,i,j,_ind,qhi,itt],
2169    iss : setify1([]),
2170    qhi : charsets_setdifference ( qh, setify1(qs)),
2171    if _is #  setify1([]) then (
2172       for i in _is do (
2173          itt : setify1(endcons( i, qs )),
2174          _ind : 0,
2175          if 0 < length ( qhi) then (
2176             for j in qhi while _ind = 0 do (
2177                if charsets_subset ( j, itt) then ( _ind : 1))),
2178          if _ind = 0 then ( iss : setify1(cons(itt, iss))))),
2179    iss
2180 )$
2181 /* subroutine for trisersub */
2182
2183 charsets_adjoina ( _is,qs,qh):=block( [ iss,i,j,_ind,qhi,itt],
2184    iss :  setify1([]),
2185    qhi : charsets_setdifference ( qh,  setify1([ qs ])),
2186    if _is #  setify1([]) then (
2187       for i in _is do (
2188          itt :  charsets_adjoin2(i, setify1(qs)),
2189          _ind : 0,
2190          if 0 < length ( qhi) then (
2191             for j in qhi while _ind = 0 do (
2192               if charsets_subset (setify1([j]), itt) then ( _ind : 1))),
2193          if _ind = 0 then (iss: charsets_adjoin2(cons(i,qs),iss)))),
2194    iss
2195 )$
2196
2197 /* if _ind = 0 then ( iss: charsets_adjoin2(charsets_adjoin2(i,qs),iss)))), */
2198 /* subroutine for irrcharser, qirrcharser and others */
2199
2200 charsets_nopsord(a,b):= if symbolp (a) then (false)
2201 else if symbolp (b) then (true) else
2202 if length( b) < length( a) then ( true) else ( false)$
2203
2204 /* remove some redundant ascending sets in cs */
2205 /*     irr=1 for irrcharser, irr=2 for qirrcharser, irr=-1 for trisersub */
2206 /*     and irr=0 for others */
2207
2208 charsets_contract ( cs,ord,irr):=block( [ i,j,mem,ts],
2209    mem :  setify1([]),
2210    ts :  setify1([]),
2211    if not (listp(cs) or charsets_setp(cs)) then break("Contract cs not list or set"),
2212    if length ( cs)  < 2 then cs
2213    else (
2214       for i : 1 thru length ( cs)  - 1 do
2215          if  not  member(i, mem) then (
2216             for j : i + 1 thru length ( cs)  do
2217                if  not  member(j, mem) then (
2218                   if charsets_linas(cs[j],ord,irr)
2219                         and charsets_contractsub ( cs[i], cs[j], ord) then (
2220                      ts :  charsets_adjoin2(cs[j], ts),
2221                      mem :  charsets_adjoin2( j, mem))
2222                   else (
2223                      if charsets_linas ( cs[i], ord, irr)
2224                         and charsets_contractsub ( cs[j], cs[i], ord) then (
2225                            ts :  charsets_adjoin2( cs[i], ts ))))),
2226       charsets_setdifference( setify1(cs), ts))
2227 )$
2228 /* check whether all polys in cs1 have remainders 0 wrt cs2  */
2229 /*      but none of their initials does: subroutine for contract */
2230
2231 charsets_contractsub( cs1,cs2,ord):=block( [ i,_is, flag:true],
2232    for i in cs1 do
2233       if charsets_premas( i, if listp(cs2) then cs2 else [cs2], ord) # 0 then (
2234                return ( flag:false)),
2235    if not flag then return(flag),
2236    _is : charsets_initialset1 ( cs1, ord),
2237    for i in _is do
2238       if charsets_premas( i, if listp(cs2) then cs2 else [cs2], ord) = 0 then (
2239                return ( flag:false)),
2240    flag
2241 )$
2242 /* check the irreducibility of as: subroutine for contract */
2243 charsets_has(x,y):=block(
2244    if atom(x) then (if x = y then true else false)
2245    else member(y,x)
2246 )$
2247 charsets_linas ( as,ord,irr):=block([i,j,n,m,nn,mm,flag],
2248    if irr = 1 then ( true)
2249    else (
2250       if irr = 2 then (
2251       flag:true,
2252       if 1 < length ( as) then (
2253          for nn:(n:1) while nn < length(as) and
2254              charsets_degreel(as[nn],ord) = 1 do(n:nn),
2255          if n < length ( as) then (
2256             for mm:(m:n+1) while mm <= length(as) and
2257                      charsets_degreel(as[mm],ord) = 1 do (m:mm),
2258             if m <= length ( as) then (return ( flag:false)))),
2259       flag)
2260       else (
2261          flag:true,
2262          for i in as do (
2263             if 1 < charsets_degreel ( i, ord) then ( return ( flag:false))),
2264          if not flag then return(false),
2265          if irr = 0 or length ( as)  < 2 then ( true)
2266          else (
2267             for i : 2 thru length ( as)  do (
2268                for j : 1 thru i do (
2269                   if charsets_has(charsets_initial(as[i],ord),
2270                         charsets_lvar(as[j],ord)) then (
2271                      return ( flag:false))))
2272                   )),
2273    flag)
2274 )$
2275 /* the quasi-irreducible char series of poly set ps */
2276
2277 charsets_qirrcharser ( ps,ord,medset):=block(
2278    [ qs,cs,_is,iss,n,ts,qsi,qhi,pi,factorset,ppi,qqi,csno,fset,mind],
2279    if  charsets_operatorp(ps[1], charsets_ListOrSet) then ( qhi :  setify1([ (ps) ]))
2280    else ( qhi :  setify1([ ps ])),
2281    qsi :  setify1([]),
2282    pi :  setify1([]),
2283    csno : 0,
2284    ppi :  setify1([]),
2285    qqi :  setify1([]),
2286    mind:if medset # charsets_basset and medset # charsets_wbasset then true else false,
2287    for n : 0 while qhi #  setify1([]) do (
2288       qhi : sort (qhi, charsets_nopsord),
2289       qs : qhi[1],
2290       ppi : charsets_select ( ppi, length ( qs) ),
2291       qqi :  charsets_union( qqi, ppi[2]),
2292       if n = 0 then ( ppi :  setify1([]))
2293       else ( ppi :  charsets_adjoin2( qs, ppi[1])),
2294       if length ( qs)  - 3 < length ( ord) then (
2295          if mind then (
2296             ml2(cs,fset, apply( assoc(medset,charsets_fcharsetsub_alist),
2297                [ qs, ord, [ setify1([]), listofvars( qs)]])),
2298             factorset : fset[1])
2299          else (
2300             ml2(cs,factorset, charsets_fcharseta ( qs, ord, medset)),
2301             if 1 < charsets_printlevel then (
2302                csno : csno + 1,
2303                print("Characteristic set produced in qirrcharser ",csno,length ( qhi),length ( qsi),length ( qs)),
2304                print ( cs))))
2305          else (
2306             if  not  mind and 2 < length ( listofvars ( cs[1]) )  then (
2307                ml2(cs,factorset, charsets_removecont( charsets_charseta( qs, ord, medset), ord)))
2308             else (
2309                ml2(cs,factorset, charsets_removecont( medset( qs, ord), ord)))),
2310          if 0 < charsets_class ( cs[1], ord)  then (
2311             ts : charsets_qirras ( cs, ord),
2312             if ts[2] = 0 then (
2313                if mind then (
2314                   if  medset = charsets_autored then (
2315                      cs : charsets_charseta(charsets_union(qs,cs),ord,charsets_charsetn))
2316                   else (
2317                      if  not  charsets_subset ( cs, qs)  then (
2318                         cs : charsets_charseta( charsets_union( qs, cs), ord, medset))),
2319                   if 1 < charsets_printlevel then (
2320                      csno : csno + 1,
2321                      print("Characteristic set produced in qirrcharser ",csno,length ( qhi),length ( qsi),length ( qs)),
2322                      print ( cs))),
2323                if  not  member(cs, pi) then (
2324                   pi :  charsets_adjoin2( cs, pi),
2325                   if 0 < charsets_class ( cs[1], ord)  then (
2326                      ts : charsets_qirras ( cs, ord),
2327                      if ts[2] = 0 then (
2328                         qsi :  charsets_adjoin2( cs, qsi),
2329                         if length ( cs) = length ( ord) then (
2330                            _is : charsets_factorps ( factorset))
2331                         else (
2332                            _is : charsets_union ( charsets_initialset ( cs, ord), charsets_factorps ( factorset) )),
2333                         iss : charsets_adjoin3 ( _is, qs, qqi)))
2334                   else (
2335                      iss : charsets_adjoin3 ( charsets_factorps ( factorset), qs, qqi)))
2336                else (
2337                   iss : charsets_adjoin3 ( charsets_factorps ( factorset), qs, qqi))),
2338             if ts[2] # 0 then (
2339                _is : charsets_union ( charsets_factorps ( factorset), ts[1]),
2340                if 1 < ts[2] then (
2341                   _is : charsets_union ( _is, charsets_initialset (  [ part(cs, makelist(i,i,1,ts[2] - 1)) ], ord) )),
2342                   iss : charsets_adjoin3 ( _is, qs, qqi)))
2343             else (
2344                iss : charsets_adjoin3( charsets_factorps( factorset), qs, qqi)),
2345       if 1 < length ( qhi) then ( qhi :  charsets_union(iss, rest(qhi)))
2346       else ( qhi : iss)),
2347       if qsi =  setify1([]) then ( [])
2348       else ( (sort ( charsets_contract ( qsi, ord, 2), charsets_lenord)))
2349 )$
2350 /* examine the irreducibility of as for qirrcharser */
2351 /* I think this function should be memoized */
2352 charsets_qirras ( as,ord):=block([ _ind,i,j,p,qs,fs,n,m,ja,last],
2353    qs :  [],
2354    _ind : 1,
2355    ja : 0,
2356    for i : 1 thru length ( as)  do (
2357       p : factor ( as[i]),
2358       fs:if mapatom(p) then [p]
2359       else (
2360          if "-" = op(p) then -charsets_dfactors(-p) else charsets_dfactors(p)),
2361       qs :  setify1([]),
2362       for j : 1 thru length ( fs)  do (
2363          if 0 < charsets_class ( fs[j], ord) then qs : charsets_adjoin2(fs[j],qs)),
2364       last:charsets_lvar ( p, ord),
2365       if charsets_degree ( qs[1], last)  < charsets_degree ( p, last) then (
2366          ja : 1,
2367          _ind : 0,
2368          return(false))),
2369    if _ind = 1 and 1 < length ( as) then (
2370       for nn:(n:1)
2371          while (n:nn) < length(as) and charsets_degreel(as[n:nn],ord)=1 do (0),
2372          if n < length ( as) then (
2373             for mm:(m : n + 1) while mm <= length( as) and charsets_degreel( as[m:mm], ord) = 1 do (0),
2374                if m <= length ( as) then (
2375                   print("Warning: factorization over algebraic field required for ics")))),
2376     [ qs, ja ]
2377 )$
2378 /* subroutine for `charsets/cfactor`  */
2379
2380 charsets_cfactorsub( f,as,ord):=block([ _ind,i,ff,fn,ffn,lind,secondLast,last],
2381    ff : num( f),
2382    if charsets_operatorp(ff,"^") then (
2383       _ind : map(lambda([f],charsets_newfactoras(f,as, ord)),ff))
2384    else if  charsets_operatorp(ff,"*") then (
2385       fn :  setify1([ (ff) ]),
2386       _ind : 1,
2387       for i : 1 thru length ( fn)  do
2388          if  charsets_operatorp(fn[i],"^") then (
2389             _ind : _ind * map(fn[i],charsets_newfactoras,as, ord))
2390          else (
2391             _ind : _ind * charsets_newfactoras ( fn[i], as, ord)))
2392    else ( _ind : charsets_newfactoras ( ff, as, ord)),
2393    lind:lcoeff( ff,charsets_lvar(ff, ord))
2394       /safeLcoeff(_ind, charsets_lvar(_ind,ord)),
2395    for i : length ( as)  thru 1 step  -  1 do (
2396       ml2(secondLast,fn,
2397           charsets_premB(num(lind),as[i],charsets_lvar(as[i],ord))),
2398       ml2(last,ffn,charsets_premB(denom(lind),as[i],charsets_lvar(as[i],ord))),
2399       lind : secondLast * ffn / last / fn),
2400    last:_ind * lind,
2401    if charsets_operatorp(last, [ "*", "^" ]) then ( simplify( last/denom ( f) ))
2402    else ( f)
2403 )$
2404 charsets_cfactor_switch : false $
2405 /* factorize poly f over algebraic number field with minimal polys in as */
2406 /*       -- a new method of Wang   */
2407 /*                          first factorize pf over smaller fields  */
2408
2409 charsets_newfactoras ( f,as,ord):=block( [ pf,pas,aas,con,vf,va,i,fn,last],
2410    vf : charsets_lvar ( f, ord),
2411    if charsets_class ( vf, ord)  <= charsets_class ( as[length ( as) ], ord)
2412           then ( return ( f))
2413    else if charsets_degree ( f, vf)  = 1 then ( return ( f)),
2414    aas :  [],
2415    con : setify1([2]),
2416    for i : 1 thru length ( as)  do (
2417       last:charsets_degreel ( as[i], ord),
2418       if 1 < last then (
2419          con : charsets_adjoin2(last,con),
2420          aas :  endcons(expand(as[i]),aas))),
2421    if aas =  [] then ( return ( f)) ,
2422    va : charsets_lvar ( aas[1], ord),
2423    if length ( aas)  = 1 and charsets_degree ( f, va)  = 0 then (
2424       last:charsets_trivial ( f, aas[1], vf, va, fn),
2425       if  charsets_operatorp(last, setify1([ "*", "^" ])) then (
2426          return ( map(last,charsets_newfactoras,as, ord) / fn))) ,
2427    pas : charsets_fcnorm ( aas, ord, 1),
2428    last:charsets_fcnormal(endcons(f,pas), ord, 1),
2429    if 1 < length (last) then ( pf : f) else ( pf : last[length(last)]),
2430    if charsets_cfactor_switch and 1 < length ( pas) then (
2431       fn : charsets_cfactor ( pf,  [ pas[length ( pas) ] ], ord),
2432       if apply(max,map([maplist(lambda([x],x),num(fn))],degree,vf))
2433          < degree ( pf, vf) then (
2434              return( charsets_cfactorsub( num( fn), pas, ord) / denom ( fn)))
2435       else (
2436          fn :charsets_cfactor( pf,makelist(pas[i],i,1,length (pas)-1), ord),
2437          if apply(max , map(lambda([x],degree(x,vf)),num( fn)))
2438             < degree ( pf, vf) then (
2439             return( charsets_cfactorsub( num( fn), pas, ord) / denom( fn))))),
2440    if  charsets_adjoin2( charsets_degreel ( pf, ord),con ) = setify1([ 2 ])
2441          and length ( pas) < 3 then (
2442             charsets_factoras ( pf, pas, ord))
2443    else (
2444       if 1 < charsets_printlevel then (
2445                print("newfactoras: factorization over algebraic field: degree",
2446                charsets_degreel ( pf, ord),terms,length ( pf))) ,
2447       charsets_con : true,
2448       charsets_newfactorassub ( pf, pas, ord))
2449 )$
2450 /* test for a trivial case --- can it be extended?       */
2451 /*                            a trivial case */
2452 /* Check the usage of last */
2453 charsets_trivial ( ff,aa,vf,va,fn):=block(
2454    [ f,a,da,df,ss,i,last,secondLast],
2455       f : expand ( ff),
2456       a : expand ( aa),
2457       da : degree ( a, va),
2458       df : degree ( f, vf),
2459       if num( f / lcoeff ( f) )  = subst(num( a / lcoeff ( a) ), ( va = vf)) and 2
2460           < length ( f)
2461           then (
2462             fn : 1,
2463             ss :  [ coeff(f,vf,da) ],
2464             for i : df - 1 thru 0 step  -  1 do
2465                ss :  [ ss[1] * va + coeff(f,vf,i), (ss) ],
2466             last:sum ( ss[i + 1] * vf ^ ( ( i - 1)), i, 1, df),
2467             return ( ( ( vf - va)) * last)),
2468       secondLast:sum ( charsets__z ^ ( ( da - i)) * va ^ i * coeff(a,va,i), i, 0, da),
2469       last:sum ( charsets__z ^ ( ( df - i)) * vf ^ i * coeff(f,vf,i), i, 0, df),
2470       last:charsets_premB ( last, secondLast, charsets__z, fn),
2471       if charsets_degree ( last, charsets__z)  = 0 then ( factor ( last))
2472          else ( f)
2473 )$
2474 /* modified pseudo-divison with multiplied initial factor as fn */
2475
2476 charsets_premB ( uu,vv,x):=block( [ r,v,dr,dv,l,t,lu,lv,gn,g],
2477    if integerp(vv / x) then [ subst(x = 0,uu),1]
2478    else (
2479       gn : 1,
2480       r : expand ( uu),
2481       dr : degree ( r, x),
2482       v : expand ( vv),
2483       dv : degree ( v, x),
2484       if dv <= dr then (
2485          l : coeff(v,x,dv),
2486          v : expand ( v - l * x ^ dv))
2487       else ( l : 1),
2488       if 1 < charsets_printlevel and 500 < length ( r)  + length ( v) then (
2489          print("pseudo-division:",[length(r),x,dr],[length(v),x,dv])),
2490          while (dv <= dr and r # 0) do (
2491             g:gcd ( lu:l, lv:coeff(r,x,dr)), lu:lu/g, lv:lv/g,
2492             t : expand ( x ^ ( ( dr - dv)) * v * lv),
2493             if dr = 0 then ( r : 0) else ( r : subst(( x ^ dr = 0),r)),
2494             r : expand ( lu * r)  - t,
2495             gn : gn * lu,
2496             dr : degree ( r, x)
2497          ),
2498       [r,gn])
2499 )$
2500 /* main subroutine for newfactoras    */
2501 /* considerably modified and improved in December 1991 */
2502
2503 charsets_newfactorassub ( f,as,ord,[args]):=block(
2504    [nord,mord,cs,ccs,ccs1,cr,i,con,ff,fff,nas,vf,die,der,ci,fs,m,n,inda,
2505     indb,_is,_CS,fmedset, ncs,cc,bb,last,flag],
2506    /* args is being used as a flag, but there are always 3 or more args */
2507    if length(args) > 0 then args:true else args:false,
2508    nas : length ( as),
2509    vf : charsets_lvar ( f, ord),
2510    if 1 < nas then (
2511       for i : 2 thru nas do (
2512          as[i]:charsets_newfactorassub(as[i],makelist(as[i],i,1,i-1),ord))),
2513    last:makelist(charsets_lvar( as[i], ord), i , 1 , nas),
2514    nord : cons( vf, last ),
2515    mord : endcons(vf, last),
2516    con : 0,
2517    for i : 2 thru length ( nord)  do
2518       con : con + charsets_degree ( f, nord[i]),
2519    if con = 0 then (
2520       if nas = 1 then (
2521          m : charsets_degree ( as[1], nord[2]),
2522          n : charsets_degree ( f, vf),
2523          if gcd ( m, n)  = 1 then ( return ( f)))),
2524    last:charsets_setdifference( listofvars(charsets_adjoin2( f, as)), setify1([nord])),
2525    if length ( last)  = 1 and args = false then (
2526       return ( charsets_tefactor ( f, as, mord,  [ last ] ) )) ,
2527    indb : true,
2528    ci :  [],
2529    cr :  [],
2530    fff : f,
2531    flag:false,
2532    while(true) do (
2533       if indb then (
2534          if con # 0 and charsets_con then (
2535             der : 0,
2536             ff : f,
2537             con : 0)
2538          else (
2539             if charsets_das #  [ false ] then (
2540                die : charsets_das[1],
2541                charsets_das :  rest(charsets_das))
2542             else ( die : charsets_die ()),
2543             der : sum ( ( ( i - die - 2)) * nord[i], i, 2, length ( nord) ),
2544             ff : expand ( subst(vf = der + vf,f))),
2545          charsets_con : false,
2546          charsets_with :  setify1([]),
2547          if 1 < charsets_printlevel then (
2548             print("Characteristic set computation:",
2549                    charsets_index(endcons(ff, as), nord))),
2550          fmedset : charsets_fqcharsetn,
2551          _is : nord[length ( nord) ],
2552          if length(as) = 1 and degree( as[1], _is) > 5 and degree( ff, _is) > 5
2553             and length( charsets_union( listofvars( as), listofvars( ff))) = 2 then (
2554             ccs : factor ( resultant ( ff, as[1], _is) ),
2555             if "+" = inpart(ccs, 0) then (return ( flag:true))
2556             else(
2557                if degree ( as[1], _is)  < degree ( ff, _is)  then (
2558                   ccs :  [ ccs, as[1] ])
2559                else ( ccs :  [ ccs, ff ])),
2560             fs :  [  setify1([]),  setify1([])  ])
2561          else (
2562             ml2(ccs,fs,fmedset(endcons(ff,as),nord,[setify1([]),setify1([])]))),
2563          _is : charsets_initialset ( ccs, nord),
2564          fs : charsets_factorps ( charsets_movefactorps ( fs[1], _is, nord) ),
2565          cs : setify1(charsets_qfactor( factor( charsets_movefactorps(
2566              ccs[1], charsets_union ( _is, fs), ord) ), nord))),
2567       if not indb or cs =  setify1([]) then (
2568          charsets_with :  setify1([]),
2569          if fs =  setify1([]) then (
2570             bb : _is[1],
2571             _is : charsets_setdifference ( _is,  setify1([ _is[1] ]) ))
2572          else (
2573             bb : fs[1],
2574             fs : charsets_setdifference ( fs,  setify1([ fs[1] ]) )),
2575             ccs : charsets_charsetn(  [ (as), ff, bb ], nord),
2576             _is : charsets_initialset ( ccs, nord),
2577             cs : setify1(charsets_qfactor ( factor ( charsets_movefactorps (
2578                 ccs[1], charsets_union ( _is, fs), ord) ), nord)),
2579             indb : true) ,
2580       if cs #  setify1([]) then (
2581          if charsets_checkwith(charsets_with,charsets_union( _is, fs)) then ( inda:true)
2582          else ( inda : false),
2583          if charsets_linearas ( ccs, nord)  or 1 < length ( cs) then (
2584             if length ( cs)  = 1 then (
2585                if inda then (
2586                   ccs1 :  [ cs[1], part(ccs, makelist(i,i,2,length ( ccs))) ])
2587                else (
2588                   ccs1 : charsets_charseta( [ (as), ff, cs[1],
2589                      part(ccs, makelist(i,i,2,length ( ccs))) ],
2590                      nord, charsets_charsetn)),
2591                while (charsets_class ( ccs1[1], nord)  = 0) do (
2592                   if fs =  setify1([]) then (
2593                      bb : _is[1],
2594                      _is : charsets_setdifference ( _is,  setify1([ _is[1] ]) ))
2595                   else (
2596                      bb : fs[1],
2597                      fs : charsets_setdifference ( fs,  setify1([ fs[1] ]) )),
2598                      ccs1 : charsets_charseta (  [ (as), ff, bb ], nord,
2599                         charsets_charsetn)),
2600                if ccs1 # ccs then (
2601                   ccs : ccs1,
2602                   cs :charsets_qfactor(factor (ccs[1]),nord))),
2603             if length ( cs)  = 1 then (
2604                if charsets_linearas ( ccs, nord) then (
2605                   cc :charsets_algcd(f,subst(vf = vf - der,ccs[1]),as,mord),
2606                   if 0 < charsets_degree ( cc, vf) then (
2607                      cs :  [ charsets_arrange ( charsets_union ( _is, fs), vf)  ],
2608                      if 1 < charsets_printlevel then (
2609                         print("A non-trivial factor found:",cc)),
2610                      if length ( cs)  = 0 then ( ci :  [ fff ], fff : 1)
2611                      else (
2612                         fff : charsets_divide ( fff, cc, as, mord),
2613                         ci :  [ cc ]),
2614                         _CS : subst(vf = vf - der,cs),
2615                         flag:false,
2616                         for i : 1 thru length ( _CS)  do (
2617                            if 2 < charsets_printlevel then (
2618                               print("Algebraic GCD:",length ( _CS),i)),
2619                            cc : charsets_algcd( fff, expand( _CS[i]), as, mord),
2620                            if 0 < charsets_degree ( cc, vf) then (
2621                               fff : charsets_divide ( fff, cc, as, mord),
2622                               if 1 < charsets_printlevel then (
2623                                  print("A non-trivial factor found:",cc)),
2624                               if charsets_degree ( cc, vf)  = 1 then (
2625                                  ci : endcons(cc, ci))
2626                               else ( cr :  endcons(cc, cr)))),
2627                         return(flag:true))
2628                   else ( indb : false)))
2629             else (
2630                ncs : length ( cs),
2631                cs : append(charsets_arrange(cs,vf),
2632                      charsets_arrange(charsets_union(_is,fs),vf)),
2633                _CS : subst(vf = vf - der,cs),
2634                for i : 1 thru length ( _CS)  do (
2635                   if 2 < charsets_printlevel then (
2636                      print("Algebraic GCD:",length ( _CS),i)) ,
2637                   cc : charsets_algcd ( fff, expand ( _CS[i]), as, mord),
2638                   if 0 < charsets_degree ( cc, vf) then (
2639                      fff : charsets_divide ( fff, cc, as, mord),
2640                      if 1 < charsets_printlevel then (
2641                         print("A non-trivial factor found:",cc)) ,
2642                      if charsets_degree(cc, vf) = 1 then (ci : endcons(cc,ci))
2643                      else (
2644                         if i <= ncs then (
2645                            con : part(ccs, makelist(i,i,2,length ( ccs))),
2646                            if inda and not charsets_vanish ( cs[i], _is) then (
2647                               con : [ cs[i], con ])
2648                            else (
2649                               con : charsets_charseta (  setify1([ (as),
2650                                  con, ff, cs[i] ]), nord, charsets_charsetn)),
2651                            if con[1] = cs[i] and charsets_linearas ( con, nord)
2652                               then ( ci : endcons( cc,ci))
2653                            else ( cr :  endcons( cc,cr)))
2654                         else ( cr :  endcons( cc,cr))))),
2655                   if 0 < length( ci) or 1 < length( cr) or length( cr)  = 1 and
2656                       charsets_degree(num( cr[1]), vf) < charsets_degree( f, vf)
2657                          then ( return(false)))))),
2658      if flag then return(f),
2659      last:charsets_degree ( fff, vf),
2660      if 1 < last then ( cr : endcons(fff, cr))
2661      else if last = 1 then ( ci : endcons(fff, ci)),
2662      if length ( ci)  = 0 then (
2663         charsets_prod(
2664            makelist(charsets_newfactorassub(cr[i],as,ord),i,1,length(cr))))
2665      else(
2666         if length ( cr)  = 0 then ( product ( ci[i], i, 1, length ( ci) ))
2667         else (
2668             product( ci[i],i,1,length(ci)) * charsets_prod(
2669                makelist(charsets_newfactorassub(cr[i],as,ord),i,1,length(cr)))))
2670 )$
2671 /* compute the GCD of f and g over the algebraic field having  */
2672 /* adjoining asc set as */
2673
2674 charsets_algcd ( f,g,as,mord):=block( [ nas,fs,last,temp],
2675    nas : length ( as)  + 1,
2676    if charsets_degree(f,mord[nas]) = 0 or charsets_degree( g, mord[nas]) = 0
2677           then ( return ( 1)) ,
2678    if 1 < charsets_printlevel then (
2679       print("GCD computation over algebraic field:",charsets_index([f,g],mord),
2680            (charsets_index ( as, mord)))),
2681    ml2(temp,fs,charsets_fcharsetnA(as,[f,g], mord, [setify1([]),setify1([])])),
2682    last:part([charsets_fcnormal(temp ,mord,2)],1),
2683    if length ( last) = nas then ( last[nas]) else ( 1)
2684 )$
2685 /* subroutine for algcd */
2686
2687 charsets_fcharsetnA ( as,ps,ord,fset1):=block([cs,rs,fset2,fset3,nas],
2688    nas : length ( as)  + 1,
2689    cs : append(as, charsets_basset ( ps, ord)),
2690    fset2 : [ fset1[1], charsets_union ( fset1[2], charsets_initialset1( cs, ord) )],
2691    if 0 < charsets_degree ( cs[nas], ord[nas]) then (
2692       rs : charsets_remseta ( charsets_minus ( ps, cs), cs, ord),
2693       /* print("before removefactor3"), */
2694       ml2(rs,fset3,charsets_removefactor ( rs, ord, fset2)))
2695    else ( return (  [ 1,fset2 ] )),
2696    if rs =  [] then return([cs, fset3])
2697    else ( charsets_fcharsetnA ( as, endcons( rs, cs[nas] ), ord, fset3))
2698 )$
2699
2700 /* division over an algebraic field with adjoining ascending set as */
2701
2702 charsets_divide ( ff,f,as,ord):=block(
2703    [ m,q],
2704       sprem ( ff, f, charsets_lvar ( ff, ord), m, q),
2705       charsets_premas ( q, as, ord)
2706 )$
2707 /* check if an ascending set is quasilinear                                     */
2708 /*                               */
2709
2710 charsets_linearas ( cs,ord):=block(
2711    [ i],
2712       if length ( cs)  = 1 then ( true)
2713          else (
2714             for i : 2 thru length ( cs)  do
2715                if 1 < charsets_degreel ( cs[i], ord) then ( return ( false)) ,
2716             true)
2717 )$
2718 /* order a set ps of polys according their degrees in x */
2719
2720 charsets_arrange ( ps,x):=block(
2721       charsets_reorderb( ps, charsets_arrangesub, x)
2722 )$
2723 /* subroutine for arrange */
2724
2725 charsets_arrangesub ( f,g,x):=block(
2726    if charsets_degree ( f, x)  < charsets_degree ( g, x) then true else false
2727 )$
2728 /* random generator for linear transformation */
2729
2730 charsets_die : random(5)+3 $
2731 /* remove polys in ps as factors from f   */
2732
2733 charsets_movefactorps ( f,ps,ord):=block( [ p,ff,i],
2734    if  not   charsets_operatorp(f, charsets_ListOrSet) then (
2735       ff : f,
2736       for p in ps do ff : first(charsets_movefactor ( ff, p, ord)),
2737       ff)
2738    else ( setify1(makelist(charsets_movefactorps(f[i],ps,ord),i,1,length(f))))
2739 )$
2740 /* check if q vansihes one poly in ps */
2741
2742 charsets_vanish ( q,ps):=block( [ p,flag:false],
2743    for p in ps do if 0 = divide ( q, p) then ( return ( flag:true)) ,
2744    flag
2745 )$
2746 /* sequence of non-constant factors of f */
2747
2748 charsets_qfactor ( f,ord):=block( [ i,last],
2749    if charsets_class ( f, ord)  = 0 then ( ( setify1([])))
2750    else if  charsets_operatorp(f,"^") then (
2751       last:part( f, 1),
2752       num( last / lcoeff ( last,ord) ))
2753    else if  charsets_operatorp(f,"*") then (
2754       maplist(lambda([p],charsets_qfactor( p, ord)), f))
2755    else ( num( f / safeLcoeff(f,ord)))
2756 )$
2757 /* sequence of non-constant (multiple) factors of f */
2758 /* there is a problem with the last line.  If qqfactor is called with
2759    a single factor, then this will not return a list.  If this function
2760    fails, then the if "*" clause needs to be replaced by a double loop,
2761    which will flatten the list returned by qqfactor. */
2762 charsets_qqfactor ( nf,ord):=block( [temp, i,last,f:nf/numfactor(nf)],
2763    if charsets_class ( f, ord)  = 0 then []
2764    else if  charsets_operatorp(f,"^") then (
2765       last:part( f, 1),
2766       makelist(num( last[i] / lcoeff ( last[i],ord) ), i, 1, part( f, 2)))
2767    else if  charsets_operatorp(f,"*") then (
2768       temp:[],
2769       for i:1 thru length(f) do (
2770          last:charsets_qqfactor( part( f, i), ord),
2771          if last # [] then temp:endcons(last,temp)),
2772       temp)
2773    else (num( f / lcoeff ( f,ord) ))
2774 )$
2775 /*  the following routines implement a heuristic procedure for poly factorization */
2776 /*  over algebraic function fields by integer substitution and solving systems */
2777 /*  of linear equations     */
2778 /* the main routine                   */
2779
2780 charsets_tefactor ( f,as,ord,var):=block(
2781    [ i,j,k,vf,nv,df,inf,ff,gg,ja,js,fs,ffs,hs,gs,sol,ci,dvar,tvar,mm,tt,yvar,das,last,secondLast],
2782       nv : length ( var),
2783       vf : ord[length ( ord) ],
2784       df : charsets_degree ( f, vf),
2785       if nv = 1 and length ( as)  = 1 then (
2786             last:charsets_prem ( f, as[1], var[1]),
2787             if last # f then (
2788                   last:factor ( last), secondLast:last,
2789                    [ last:charsets_qqfactor ( last,  [ vf ] )  ],
2790                   if  charsets_operatorp(secondLast, setify1([ "*", "^" ])) and max ( last:map(charsets_degree ( last[i], vf)))  < df then (
2791                         last:map(part( part(  [ charsets_fcnormal (  [ as[1], last[i] ], ord, 2)  ], 1), 2)),
2792                         return ( charsets_prod ( map(last,charsets_newfactoras,as, ord)) )))),
2793       inf : lcoeff ( f, vf), secondLast:last,last:inf,
2794       for i : 1 thru nv do charsets_m[i] : 1,
2795       js :  [],
2796       fs :  [],
2797       charsets_dasA :  [ 1, -1, 2, -2, 3, -3, false ],
2798       charsets_dieA : random(20*nv)-10*nv,
2799       das : lambda ([], block(
2800          if charsets_dasA #  [ false ] then (
2801             charsets_dasA[1],
2802             charsets_dasA :  [ part(charsets_dasA, makelist(i,i,2,length ( charsets_dasA))) ],
2803             secondLast)
2804          else ( charsets_dieA ()))),
2805       tvar : charsets_noterms ( nv, charsets_degree ( f, var) ),
2806       dvar : 0,
2807       ci :  setify1([]),
2808       gg : _y1,
2809       yvar :  setify1([ _y1 ]),
2810       for mm : 1 while ci =  setify1([])  and dvar < tvar do (
2811          dvar : charsets_noterms ( nv, mm),
2812          while (length ( js)  < dvar) do (
2813             /* I am assuming that the following is a assignment */
2814             /* sol : map(var[i] = das ()), */
2815             for x in var do x:das(), sol:var,
2816             if subst(inf, sol) # 0 and charsets_isirr ( subst(as, sol), ord)
2817                 then (
2818                   js :  [ (js),  setify1([ sol ])  ],
2819 /* check this line carefully */
2820                   ff : setify1(charsets_qfactor ( charsets_cfactor ( subst(f, ( sol)), subst(as, sol), ord),  [ vf ] )),
2821                   if max ( map(charsets_degree ( ff[i], vf))) = df
2822                       then (
2823                         return ( f))
2824                      else (
2825                         fs :  [ (fs), ff ]))
2826         ),
2827          apply("+",var),
2828          coeffs ( expand ( % ^ mm), var, tt),
2829          tt :  [ tt ],
2830          length ( gg),
2831          /* where is _y and _y1 defined */
2832          gg : gg + sum ( _y [ % + i] * tt[i], i, 1, length ( tt) ),
2833          yvar : charsets_union ( yvar,  setify1([ map(_y [ %% + i]) ]) ),
2834          hs :  setify1([]),
2835          ffs : fs,
2836          for j : 1 thru length ( fs[1])  do (
2837             gs :  [ fs[1][j] ],
2838             if 1 < length ( fs)
2839                 then (
2840                   for i : 2 thru length ( fs)  do (
2841                      charsets_getclose ( fs[1][j],  [ (fs[i]) ], ord, ja),
2842                      if % = FAIR
2843                          then (
2844                            if 1 < charsets_printlevel
2845                                then (
2846                                  print("Heuristic tefactor failed")),
2847                            return ( charsets_newfactorassub ( f, as, ord, 0) ))
2848                         else (
2849                            gs :  [ (gs), % ]),
2850                      charsets_setdifference ( fs[i],  setify1([ fs[i][ja] ]) ),
2851                      if i = length ( fs)
2852                          then (
2853                            fs :  [ part(fs, makelist(i,i,1,i - 1)), % ])
2854                         else (
2855                            fs :  [ part(fs, makelist(i,i,1,i - 1)), %, part(fs, makelist(i,i,i+ 1,length ( fs))) ]))),
2856             sol :  setify1([ map(subst(gg, ( (js[k]))) - gs[k]) ]),
2857             sol :  setify1([ solve ( sol, yvar)  ]),
2858             if sol #  setify1([])
2859                 then (
2860                   hs : charsets_union ( hs,  setify1([ expand ( subst(gg, ( (sol))))  ]) ))),
2861          fs : ffs,
2862          ff : f,
2863          for j in hs do (
2864             charsets_divideA ( ff, j, as, ord),
2865             if % # false
2866                 then (
2867                   ff : %,
2868                   ci : charsets_union ( ci,  setify1([ j ]) ),
2869                   if 1 < charsets_printlevel
2870                       then (
2871                         print("A factor found:",j)))),
2872          if ci =  setify1([])  and mm <= 1 and _help # true and length ( ffs[1])  ^ length ( ffs)
2873              <= 128
2874              then (
2875                gs : charsets_getall ( ffs),
2876                for i : 1 thru length ( gs) while length ( ci)  <= length ( fs[1]) do (
2877                   sol :  setify1([ map(subst(gg, ( (js[k]))) - gs[i][k]) ]),
2878                   sol :  setify1([ solve ( sol, yvar)  ]),
2879                   if 1 < charsets_printlevel and sol =  setify1([])
2880                       then (
2881                         print(sol,yvar)),
2882                   if sol #  setify1([])
2883                       then (
2884                         sol : expand ( subst(gg, ( (sol)))),
2885                         charsets_divideA ( ff, sol, as, ord),
2886                         if % # false
2887                             then (
2888                               ff : %,
2889                               ci : charsets_union ( ci,  setify1([ sol ]) ),
2890                               if 1 < charsets_printlevel
2891                                   then (
2892                                     print("A non-trivial factor found:",j))))))),
2893       if ci #  setify1([]) then (
2894             ci : charsets_union ( ci,  setify1([ ff ]) ),
2895             charsets_prod ( map(ci,charsets_newfactoras,as, ord)))
2896          else (
2897             if 1 < charsets_printlevel
2898                 then (
2899                   print("Heuristic tefactor failed")),
2900             charsets_newfactorassub ( f, as, ord, 0))
2901 )$
2902 /* numbers of maximal terms in a poly of total degree d in n variables */
2903
2904 charsets_noterms ( n,d):=block([ i,j],
2905    sum ( product ( n + j - 1, j, 1, i)  / i!, i, 1, d)  + 1
2906 )$
2907 /* check if an ascending set as is irreducible */
2908
2909 charsets_isirr ( as,ord):=block( [ xa,fs,f,as1],
2910    if length ( as)  = 1 then (
2911       xa : charsets_lvar ( as[1], ord),
2912       if xa = 0 then ( false)
2913       else if charsets_degree ( as[1], xa)  = 1 then ( true)
2914          else (
2915             fs :  setify1(charsets_qfactor ( factor ( as[1]),  [ xa ] )),
2916             if length ( fs)  = 1 then ( true) else ( false)))
2917       else (
2918          as1 :  [ part(as, makelist(i,i,1,length ( as)  - 1)) ],
2919          if  not  charsets_isirr ( as1, ord) then ( false)
2920          else (
2921             f : as[length ( as) ],
2922             xa : charsets_lvar ( f, ord),
2923             if xa = 0 then ( false)
2924             else if charsets_degree ( f, xa)  = 1 then ( true)
2925             else (
2926                fs : setify1(charsets_qfactor(charsets_cfactor(f,as1,ord),[xa])),
2927                if length ( fs)  = 1 then ( true) else ( false))))
2928 )$
2929 /* get all possible combinations (used for tefactor) */
2930
2931 charsets_getall ( fs):=block(
2932    [ gs,i,j,nf],
2933       if length ( fs)  = 1
2934           then (
2935              setify1([ map( [ fs[1][i] ] /* $ ( ( i = 1 .. length ( fs[1]) ))*/) ]))
2936
2937          else (
2938             nf : length ( fs),
2939             gs :  setify1([ part(fs, makelist(i,i,1,nf - 1)) ]),
2940             gs : charsets_getall ( gs),
2941              setify1([ map(map( [ (gs[i]), fs[nf][j] ] /* $ ( ( j = 1 .. length ( fs[nf]) ))*/)/* $ ( ( i
2942                 = 1 .. length ( gs) ))*/) ]))
2943 )$
2944 /* select a poly in fs closest to g */
2945
2946 charsets_getclose ( g,fs,var,ja):=block(
2947    [ i,j,gs,hs,dg,nv,vv,jaa,jbb,ts,cv,rr,chs,df,last],
2948    if charsets_class ( g, var)  = 0 then ( return ( FAIR)),
2949    if length ( fs)  = 1 then ( ja : 1, return ( fs[1])),
2950    if length ( var)  = 1 and _help # true then (
2951       gs :  [],
2952       jaa :  [],
2953       dg : charsets_degree ( g, var[1]),
2954       for i : 1 thru length ( fs)  do (
2955          if fs[i] = g or fs[i] =  -  g then ( ja : i, return ( fs[i]))
2956          else if charsets_degree ( fs[i], var[1])  = dg then (
2957             jaa :  [ (jaa), i ],
2958             gs :  [ (gs), fs[i] ])
2959      ),
2960       if length ( jaa)  = 1 then ( ja : jaa[1], return ( gs[1])),
2961       hs :  [],
2962       jbb :  [],
2963       cv :  [ coeffs ( expand ( g), var[1], dg)  ],
2964       for i : 1 thru length ( jaa)  do (
2965          coeffs ( expand ( gs[i]), var[1], df),
2966          if  setify1([ df ])  =  setify1([ dg ])  then (
2967             jbb :  [ (jbb), jaa[i] ],
2968             hs :  [ (hs), gs[i] ])),
2969          if length ( jbb)  = 1 then ( ja : jbb[1], return ( hs[1])),
2970          gs :  [],
2971          jaa :  [],
2972          dg :  [ map(sign, cv) ],
2973          for i : 1 thru length ( jbb)  do (
2974             ts :  [ coeff(expand ( hs[i]),[var[1]])],
2975             ts :  [ map(sign, ts) ],
2976             if charsets_close ( ts, dg)  = 0 then (
2977                jaa :  [ (jaa), jbb[i] ],
2978                gs :  [ (gs), hs[i] ])),
2979          if length ( gs)  = 1 then ( ja : jaa[1], return ( gs[1])),
2980          gs :  [],
2981          jaa :  [],
2982          for i : 1 thru length ( jbb)  do (
2983             ts :  [ coeff(expand ( hs[i]),[var[1]])],
2984             ts :  [ map(sign ( ts[j])) ],
2985             if charsets_close ( ts, dg)  = 1 then (
2986                jaa :  [ (jaa), jbb[i] ],
2987                gs :  [ (gs), hs[i] ])),
2988          if length ( gs)  = 1 then (ja : jaa[1], return ( gs[1])),
2989          gs :  [],
2990          jaa :  [],
2991          for i : 1 thru length ( jbb)  do (
2992             ts :  [ coeff(expand ( hs[i]),[var[1]])],
2993             ts :  [ map(sign, ts) ],
2994             if charsets_close ( ts, dg)  = 2 then (
2995                jaa :  [ (jaa), jbb[i] ],
2996                gs :  [ (gs), hs[i] ])),
2997          if length ( gs)  = 1 then ( ja : jaa[1], return ( gs[1]))
2998          else ( return ( FAIL)))
2999       else (
3000          if _help # true then (
3001             nv : length ( var),
3002             vv : var[nv],
3003             gs :  [],
3004             jaa :  [],
3005             dg : charsets_degree ( g, vv),
3006             for i : 1 thru length ( fs)  do(
3007                if fs[i] = g or fs[i] = -g then ( ja : i, return ( fs[i]))
3008                else if charsets_degree ( fs[i], vv)  = dg then (
3009                   jaa :  [ (jaa), i ],
3010                   gs :  [ (gs), fs[i] ])),
3011             if length ( jaa)  = 1 then ( ja : jaa[1], return ( gs[1])),
3012             hs :  [],
3013             jbb :  [],
3014             hs :  [],
3015             chs :  [],
3016             /* I think here is the calls to coeff with 3 args */
3017             cv :  [ coeffs ( expand ( g), vv, dg)  ],
3018             for i : 1 thru length ( gs)  do(
3019                chs :  [ (chs),  [ coeffs ( expand ( gs[i]), vv, df)  ]  ],
3020                if  setify1([ df ])  =  setify1([ dg ]) then (
3021                   jbb :  [ (jbb), jaa[i] ],
3022                   hs :  [ (hs), gs[i] ])),
3023             if length ( jbb)  = 1 then ( ja : jbb[1], return ( hs[1])),
3024             for i : 1 thru length (  [ dg ] )  do(
3025                ts :  [ map(chs[j][i]/* $ ( ( j = 1 .. length ( chs) ))*/) ],
3026                charsets_getclose ( cv[i], ts,  [ map(var[j]/* $ ( ( j = 1 .. nv - 1))*/) ], jbb),
3027                if % = FAIR then ( return ( FAIR))
3028                else if % # FAIL then ( ja : jaa[jbb], return ( hs[jbb]))))
3029          else (
3030             hs : fs,
3031             jaa : makelist(i, i, 1,length ( hs) )),
3032       if hs #  []  then (
3033             return ( FAIR),
3034             print("Please help to choose one polynomial in the list",hs),
3035             print("which is closest to the polynomial",g),
3036             rr : readstat ("and enter the polynomial number in the list:" ),
3037             ja : jaa[rr],
3038             return ( hs[rr]))
3039       else ( ja : 1, return ( fs[1])))
3040 )$
3041 /* a subroutine for getclose */
3042
3043 charsets_close ( ps,qs):=block(
3044    [ i,m],
3045       if ps = qs
3046           then (
3047             0)
3048
3049          else (
3050             m : 0,
3051             for i : 1 thru length ( ps)  do
3052                if ps[i] # qs[i]
3053                    then (
3054                      m : m + 1)
3055 ,
3056             m)
3057 )$
3058 /* division over an algebraic field with adjoining ascending set as */
3059
3060 charsets_divideA ( ff,f,as,ord):=block(
3061    [ m,q,last],
3062       if charsets_class ( ff, ord)  # charsets_class ( f, ord)
3063           then (
3064             return ( false))
3065 ,
3066       last:sprem ( ff, f, charsets_lvar ( ff, ord), m, q),
3067       if charsets_premas (last, as, ord)  # 0
3068           then (
3069             return ( false))
3070 ,
3071       charsets_premas ( q, as, ord)
3072 )$
3073 /* factorize poly f over algebraic number field with minimal polys in as */
3074 /*       -- Hu-Wang's method */
3075
3076 charsets_factoras ( pf,pas,ord,[args]):=block(
3077    [df,r,s,ff,i,t,fact,gg,hh,fg,z,_ind,as,nas,f,vf,sol,con,m,n,mord,nord,last],
3078    if not polynomialp(pf,ord,
3079        lambda([x],integerp(x) or (not constantp(x) and lfreeof(ord,x)))) then (
3080       return(charsets_factoras(num(pf),pas,ord)/denom(pf))),
3081    /* args is being used as a flag, but there are always 3 or more args */
3082    if length(args) > 0 then args:true else args:false,
3083    nas : length ( pas),
3084    vf : charsets_lvar ( pf, ord),
3085    if 1 < nas then (
3086       for i : 2 thru nas do (
3087          pas[i]:charsets_factoras(pas[i],makelist(pas[i],i,1,i-1),ord))),
3088
3089 /* Not quite.  Assuming that pas is a list, it assigns pas[2] as the value
3090 of `charsets/factoras`(pas[2],[pas[1]],ord).  And then for i=3,
3091 it assigns pas[3] as the value of
3092    `charsets/factoras`(pas[3],[pas[1],pas[2]],ord).
3093 These values are stored in a "remember table", so the next time the
3094 procedure is called with those arguments the value will be taken from
3095 the table rather than being recomputed.
3096 */
3097    last:makelist(charsets_lvar( pas[i], ord), i , 1 , nas),
3098    nord :  cons( vf, last ),
3099    mord :  endcons(vf, last),
3100    con : sum(charsets_degree ( pf, nord[i]),i,2,length(nord)),
3101    if con = 0 then (
3102       if nas = 1 then (
3103          m : charsets_degree ( pas[1], nord[2]),
3104          n : charsets_degree ( pf, vf),
3105          if gcd ( m, n)  = 1 then ( return ( pf)))),
3106    last:charsets_setdifference( listofvars( charsets_adjoin2( pf, pas)), setify1(nord)),
3107    if length ( last)  = 1 and args = false then (
3108       /* here it is assumed that sets are implemented as lists */
3109       return ( charsets_tsfactor ( pf, pas, mord, last) )),
3110    _ind : 0,
3111    as :  [],
3112    r : 0,
3113    f : expand ( pf),
3114    for i : 1 thru length ( pas)  do (
3115       df : charsets_degreel ( pas[i], ord),
3116       if 1 < df then (
3117          r : r + 1,
3118          charsets_m[r] : df,
3119          as :  endcons(expand ( pas[i]),as))),
3120    z : charsets_lvar ( f, ord),
3121    df : charsets_degree ( f, z),
3122    if df = 1 then ( f)
3123    else if r = 0 then ( f)
3124    else (
3125       if 1 < charsets_printlevel then (
3126          print("factoras: factorization over algebraic field -- degree ",
3127              charsets_degreel ( f, ord))),
3128       for s : 1 thru fix( (1/2) * df) while _ind = 0 do (
3129          for i : 1 thru s do block([g],
3130             last : makelist(concat(charsets_k,t),t,1,r),
3131             g:concat(charsets_g,i),
3132             g::charsets_summ( arraymake(charsets_g,cons(i,last))*product(
3133                 charsets_lvar(as[t],ord)^concat(charsets_k,t),t,1,r),r)),
3134          for i : 1 thru df - s do block([h],
3135             last : makelist(concat(charsets_k,t),t,1,r),
3136             h:concat(charsets_h,i),
3137             h:: charsets_summ( arraymake(charsets_h,cons(i,last))*product(
3138                  charsets_lvar(as[t],ord)^concat(charsets_k,t),t,1,r),r)),
3139          last:concat(charsets_g,0),
3140          last :: 1,
3141          last:concat(charsets_h,0),
3142          last :: 1,
3143          gg : sum( eval(concat(charsets_g,i))* z ^(s - i), i, 0, s),
3144          hh : sum( eval(concat(charsets_h,i))* z ^(df - s - i), i, 0, df - s),
3145          ff : f - lcoeff ( expand ( f), z)  * expand ( gg * hh),
3146          ff : expand ( charsets_premas ( ff, as, ord) ),
3147          fact :  setify1([]),
3148          for i : 0 thru df - 1 do
3149             fact:charsets_union(charsets_coeff(as,setify1([coeff(ff,z,i)]),r,ord),fact),
3150          sol :  [ charsets_solveps ( fact, charsets_getvars ( fact) )  ],
3151          if sol =  [ fair ]  then (
3152             if 1 < charsets_printlevel then (
3153                print("factoras changes to newfactoras")),
3154             charsets_das :  [ -1, 1, -2, 2, -3, false ],
3155             charsets_con : true,
3156             return ( charsets_newfactorassub( pf, pas, ord) )),
3157          fg : f,
3158          if sol #  []  then (
3159             fg:subst(sol[1],gg)*charsets_factoras(num(subst(sol[1],hh)),as,ord),
3160             _ind : 1)),
3161       num( fg))
3162 )$
3163 /* the following routine implements some heuristics for verifying the */
3164 /* irreducibilty of polynomials over algebraic function fields by  */
3165 /* integer substitution */
3166
3167 charsets_tsfactor ( f,as,ord,var):=block(
3168    [ i,vf,nv,df,_inf,ff,sol,das,secondLast,last],
3169    nv : length ( var),
3170    vf : ord[length ( ord) ],
3171    df : charsets_degree ( f, vf),
3172    if nv = 1 and length ( as)  = 1 then (
3173       last:charsets_prem ( f, as[1], var[1]),
3174       if last # f then (
3175          last:factor ( last), secondLast:last/numfactor(last),
3176          last:charsets_qqfactor( last,  [ vf ] ),
3177          if  charsets_operatorp(secondLast, [ "*", "^" ]) and
3178             apply(max, map(lambda([x],charsets_degree(x,vf)),last)) < df then (
3179                 last:makelist(
3180                    second(charsets_fcnormal([as[1],last[i]],ord,2)),
3181                    i,1,length(last)),
3182             return(charsets_prod(map(lambda([x],charsets_factoras(x,as, ord)),
3183                last)))))),
3184    _inf : lcoeff ( expand ( f), vf),
3185    /* das : rand (  -2 * nv .. 3 * nv + length ( ord)), check next line */
3186    sol : makelist(var[i] = i + 1,i, 1,nv),
3187    while (lratsubst(sol,_inf)=0 or not charsets_isirr(lratsubst(sol,as),ord)) do
3188       sol : makelist(var[i] = rand (5*nv +length(ord))-  2 * nv,i,1,nv),
3189    ff :  setify1(charsets_qfactor(
3190       charsets_cfactor( lratsubst(sol,f), lratsubst(sol,as), ord), [ vf ] )),
3191    if apply(max,map(lambda([x],charsets_degree(x,vf)),ff)) = df then return(f),
3192    if 1 < charsets_printlevel then (
3193             print("Heuristic tsfactor failed")),
3194    charsets_factoras ( f, as, ord, 0)
3195 )$
3196 /* subroutine for factoras */
3197 /* I am not sure if I translated this one correctly */
3198 charsets_summ( ss,r):=block(
3199    [temp,summ,m,_ind],
3200    if r = 1 then (
3201       m:arrayapply(charsets_m,[r])-1,
3202       temp:0,
3203       for i:0 thru m do (
3204          _ind:concat(charsets_k,r),
3205          temp:temp + subst(i,_ind,ss)),
3206       temp)
3207    else (
3208       temp:0,
3209       for i:0 thru charsets_m[r-1] do (
3210          summ:charsets_summ( ss, r - 1),
3211          temp:temp+summ[i]))
3212 )$
3213 /* subroutine for factoras */
3214
3215 charsets_coeff(as,ss,r,ord):=block([k,i,j,qs, temp],
3216    qs : ss,
3217    temp : setify1([]),
3218    for j : r thru 1 step -1 do (
3219       for i : 1 thru length(qs) do (
3220          for k : 0 thru charsets_m[j] -1 do (
3221             temp : charsets_adjoin2(coeff(qs[i],charsets_lvar( as[j], ord),k),temp)))),
3222    charsets_setdifference( temp, setify1([0]))
3223 )$
3224 /* subroutine for factoras */
3225
3226 charsets_getvars ( as):=block([ _ind,ind1,i],
3227    if charsets_operatorp(as, charsets_ListOrSet) then
3228       apply(charsets_union,map(charsets_getvars,as))
3229    else (
3230       _ind :  setify1([]),
3231       ind1 : listofvars ( as),
3232       for i in ind1 do
3233          if subvarp(i) then (
3234             if part( i, 0) = charsets_g or part( i, 0) = charsets_h then (
3235                _ind : charsets_adjoin2(i, _ind))),
3236       _ind)
3237 )$
3238 /* find rational zeros of poly set ps */
3239
3240 charsets_solveps ( ps,lst):=block(
3241    [ cs,ord,sol,j,phi,qs,qs1,n,factorset,last,temp,flag],
3242    if 1 < charsets_printlevel then
3243       print("solveps: trying rational solutions of equations:",
3244          length(ps),length(lst)),
3245    if 3 < charsets_printlevel then print(charsets_index( ps, lst)),
3246    ord : charsets_reorder(lst, charsets_degord, ps),
3247    sol :  setify1([]),
3248    ml2(cs,factorset, charsets_fqcharsetn( ps, ord, [setify1([]),setify1([])])),
3249    factorset : factorset[1],
3250    phi :  setify1([ ps ]),
3251    flag:false,
3252    for n : 1 while phi #  setify1([]) do (
3253       if 6 < n + length( phi) or sol = setify1([FAIR]) then return( flag:FAIR),
3254       if sol # setify1([])  then return (false),
3255       if 1 < n then (
3256          cs : charsets_charseta ( phi[1], ord, charsets_wcharsetn),
3257          factorset :  setify1([])),
3258       ml2(temp,qs1,charsets_solveasr( cs, ord)),
3259       sol : charsets_union(temp, sol),
3260       if n = 1 then (sol : charsets_verify ( sol, ps, ord)),
3261       qs : charsets_union ( charsets_factorps ( qs1), charsets_factorps ( factorset) ),
3262       last:charsets_qbasset ( qs, ord),
3263       qs :  charsets_union(last, charsets_setdifference ( qs,  setify1([ (last) ]) )),
3264       if qs #  setify1([])  then (
3265          if 1 < length ( phi)  then (
3266             phi :  setify1([ part(phi, makelist(i,i,2,length ( phi))), map( [ (phi[1]), qs[j] ] /* $ ( ( j
3267                          = 1 .. length ( qs) ))*/) ]))
3268
3269          else (
3270             phi :  setify1([ map( [ (phi[1]), qs[j] ] /* $ ( ( j
3271                          = 1 .. length ( qs) ))*/) ])))
3272
3273          else (
3274             if 1 < length ( phi)  then (
3275                phi :  setify1([ part(phi, makelist(i,i,2,length ( phi))) ]))
3276
3277             else ( phi :  setify1([])))),
3278    if flag # false then return(flag),
3279    if sol =  setify1([])  then (
3280       ( setify1([])))
3281    else ( (sol))
3282 )$
3283 /* subroutine for solveps */
3284
3285 charsets_verify ( sol,ps,ord):=block( [ i,j,sss,last,retval],
3286    retval:charsets_setify([]),
3287    for soli in sol  do
3288       if setify1(simplify( subst(soli,ps))) = setify1([0]) then
3289          return( retval:setify1(soli))
3290       else (
3291          last:charsets_setify([]),
3292          for solij in soli do last:charsets_adjoin2(first(solij)- second(solij),last),
3293          sss : setify1(charsets_solveps( charsets_union ( ps,last), ord)),
3294          if sss #  setify1([])  then return( retval:sss)),
3295    retval
3296 )$
3297 /* prepare a list of triangular forms from poly set ps */
3298
3299 charsets_trisersub ( ps,ord):=block(
3300    [ qs,cs,iss,n,i,qhi,qsi,factorset,csno,ppi,qqi,_ind,mem,fset4],
3301    _ind : 0,
3302    for i : 1 thru length ( ps)  do(
3303       if length ( expand ( ps[i]) )  < 3 then (
3304          _ind : 1,
3305          return(false))),
3306    if _ind = 1 then (
3307       ml2( cs, factorset, charsets_fcharseta( ps, ord, charsets_charsetn)))
3308    else (
3309       ml2( cs, factorset, charsets_fcharseta( ps, ord, charsets_qcharsetn))),
3310    qhi :  setify1([setify1(ps)]),
3311    qsi :  setify1([]),
3312    csno : 0,
3313    ppi :  setify1([]),
3314    qqi :  setify1([]),
3315    for n : 0  while qhi #  setify1([])   do (
3316       qhi : sort ( qhi, charsets_nopsord),
3317       qs : qhi[1],
3318       ppi : charsets_select ( ppi, length ( qs) ),
3319       qqi :  charsets_union(qqi,ppi[2]),
3320       if n = 0 then ( ppi :  setify1([]))
3321       else (
3322          ppi : charsets_union( setify1(qs), setify1(ppi[1])),
3323          _ind : 0,
3324          for i : 1 thru length ( ps)  do(
3325             if length ( expand ( ps[i])) < 3 then ( _ind : 1,return(false))),
3326          if _ind = 1 then (
3327             cs:charsets_nopower(charsets_charseta( qs, ord, charsets_charsetn)),
3328             factorset :  setify1([]))
3329          else ( fset4: qs,
3330            fset4: if (listp(qs) and (qs # []) and (listp(first(qs))))
3331            then first(qs) else qs,
3332            if (qs # mem) and (4 < charsets_degree (fset4[1], ord))  then (
3333              ml2(cs,factorset,charsets_fcharseta(qs,ord,charsets_qcharsetn)))
3334            else (
3335              if (length ( qs)  - 3) < length ( ord)  then (
3336                ml2(cs,factorset, charsets_fcharseta(qs,ord,charsets_wcharsetn)))
3337              else (
3338                cs : charsets_nopower (
3339                  charsets_charseta ( qs, ord, charsets_wcharsetn)),
3340                factorset :  setify1([]))))),
3341        mem : qs,
3342        if 1 < charsets_printlevel then (
3343          csno : csno + 1,
3344          print("Characteristic set produced",csno,length(qhi),
3345            length(qsi),length (qs),cs)),
3346        if 0 < charsets_class ( cs, ord)  then (
3347          iss : charsets_initialset ( cs[1], ord),
3348          if charsets_simpa ( iss, cs, ord)  # 0 then (
3349            qsi : charsets_union([cs],qsi)),
3350          iss : charsets_union ( iss, charsets_factorps ( factorset)))
3351        else ( iss : charsets_factorps ( factorset)),
3352       iss : charsets_adjoina ( iss, qs, qqi),
3353       if 1 < length ( qhi) then ( qhi : charsets_union(setify1([iss]), rest(qhi)))
3354       else ( qhi : iss)
3355    ),
3356    if qsi =  setify1([])  then ( [])
3357    else (sort ( charsets_contract ( qsi, ord,  -  1), charsets_lenord))
3358 )$
3359
3360
3361
3362 /* find zeros of ascending set as */
3363
3364 charsets_solveas ( cs,ord):=block( [ _is,ss,sol,solm,i,j,k,temp],
3365    sol:map("[",sort(solve(cs[1],charsets_lvar(cs[1],ord)),ordergreatp)),
3366    if 1 < length ( cs)  then (
3367       for i : 2 thru length ( cs)  do (
3368          _is : charsets_initial ( cs[i], ord),
3369          solm :  setify1([]),
3370          for j : 1 thru length ( sol)  do (
3371             ss : map("[",sort(solve(setify1([subst(sol[j],cs[i])]),
3372                 charsets_lvar ( cs[i], ord)),ordergreatp)),
3373             for k : 1 thru length ( ss)  do (
3374                temp:subst(ss[k],_is),
3375                temp:subst(sol[j],temp),
3376                if temp # 0 then (
3377                   solm : charsets_adjoin2(append(sol[j], ss[k]),solm)))),
3378          sol : solm)),
3379    sol
3380 )$
3381 /* find rational zeros of ascending set cs */
3382 /*
3383 fx(q):=block([c],
3384    c:numfactor(q),
3385    q:q/c,
3386    num(c)*num(q/lcoeff(q))
3387 )$
3388 */
3389 charsets_solveasr ( cs,ord):=block([ _is,ss,ts,sol,solm,i,j,k],
3390    ts :  setify1([]),
3391    if 0 < charsets_class ( cs[1], ord) then (
3392       sol :  setify1(charsets_solvel( cs[1], charsets_lvar( cs[1], ord))),
3393       if 1 < length ( cs) then (
3394          for i : 2 thru length ( cs) do
3395             if 1 <= length ( sol) then (
3396                _is : charsets_initial ( cs[i], ord),
3397                solm :  setify1([]),
3398                for j : 1 thru length ( sol) do
3399                   if simplify( subst(flatten(sol[j]),_is)) = 0 then (
3400                      ts : setify1(cons( _is, ts)))
3401                   else (
3402                      ss :setify1(charsets_solvel( subst(flatten(sol[j]),cs[i]),
3403                         charsets_lvar( cs[i], ord))),
3404                      for k in ss do
3405                         solm : charsets_adjoin2(charsets_union(k, sol[j]),solm)),
3406                      sol : solm)
3407                else ( return(false))))
3408       else ( sol :  setify1([])),
3409       [sol,ts]
3410 )$
3411
3412 /* find rational zeros of polynomial f wrt x: subroutine for solveasr */
3413
3414 charsets_solvel ( f,x):=block( [ g,i,sol],
3415    sol :  setify1([]),
3416    if length( listofvars ( f) )  = 1 then (
3417       g : charsets_getfactor( f, x),
3418       for i in g do
3419          sol :  charsets_union(solve (  setify1([ i ]),  setify1([ x ]) ),sol))
3420    else (
3421       g : charsets_factorps(  setify1([ num( f)  ]) ),
3422       for i in g do
3423          if charsets_degree ( i, x)  = 1 then (
3424             sol : charsets_union(solve( setify1([ i ]), setify1([ x ])),sol))),
3425    (sol)
3426 )$
3427 /* find a list of distinct linear factors of univariate poly f */
3428 charsets_getfactor ( ff,x):=block( [ last,f,i],
3429    f: factor(ratsimp(ff/gcd(ff,diff(ff,x)))),
3430    if charsets_degree(f,x) = 1 then [f]
3431    else (
3432       last:setify1([]),
3433       for i in f do if charsets_degree(i,x) = 1 then last:charsets_adjoin2(i,last),
3434       last)
3435 )$
3436 /*
3437 charsets_getfactor ( f,x):=block( [ q,qs,j],
3438    q : charsets_getfact ( f, x),
3439    qs : setify1([]),
3440    if charsets_operatorp(q,"*") then (
3441       for j : 1 thru length ( q)  do
3442          if not integerp(part( q, j)) then (
3443             if charsets_operatorp(part( q, j),"^") then (
3444                qs : setify1([ qs, num( part( part( q, j), 1)
3445                             / lcoeff ( part( part( q, j), 1)) )  ]))
3446             else (
3447                qs : setify1([ qs, num( part( q, j) / lcoeff( part( q, j)))]))))
3448       else if  charsets_operatorp(q,"^") then (
3449          qs :  setify1([ qs, num( part( q, 1) / lcoeff ( part( q, 1)) )]))
3450       else (
3451          if not integerp(q) then (
3452             qs : setify1([ num( q / lcoeff ( q) ), qs ]))),
3453     [ qs ]
3454 )$
3455 */
3456 /* find the product of linear factors of univar poly f */
3457
3458 charsets_getfact ( ff,x):=block([i,f,last],
3459    if charsets_degree ( ff, x)  = 1 then return ( ff),
3460    f: factor(ratsimp(ff/gcd(ff,diff(ff,x)))),
3461    last:setify1([]),
3462    for i in f do
3463       if charsets_degree(i,x) = 1 then last:charsets_adjoin2(i,last),
3464    apply("*",last)
3465 )$
3466 /*  irreducible decomposition of algebraic variety defined by ps */
3467
3468 charsets_irrvardec ( ps,ord,medset,[args]):=block(
3469    [ phi,psi,i,j,mem,qq,qs],
3470    if length(args) > 0 then args:true else args:false,
3471    qq : length ( ps),
3472    mem :  setify1([]),
3473    if 1 < charsets_printlevel then ( print("Variable order chosen:",ord)),
3474    if not args then psi : charsets_irrcharser ( ps, ord, medset)
3475    else (
3476       psi : charsets_exirrcharser (  [ ps, 1 ], ord, medset),
3477       if psi #  [  []  ] then (
3478          phi : psi,
3479          psi : ( []),
3480          for i in phi do
3481             if listp(i[1]) then
3482                psi : append(psi, i[1])
3483             else psi:endcons(i, psi),
3484          psi :  [ psi ])),
3485    phi :  [],
3486    for i : 1 thru length ( psi)  do (
3487       if length ( psi[i])  <= qq then phi :  endcons(psi[i],phi)),
3488 break("phi before sort",phi),
3489    phi : sort ( phi, lambda ([ a,b ], is( length ( a)  < length ( b)))),
3490 break("phi",phi),
3491    psi :  [],
3492    if phi #  [  []  ]  then (
3493       if length ( phi)  = 1 and args = true then (
3494          return (  [ (ps) ] )),
3495       if 1 < charsets_printlevel then (
3496          print("ivd: ics finished at",time (),length ( phi))),
3497       for i : 1 thru length ( phi)  do(
3498          if  not  member(i, mem) then (
3499             qs : charsets_primebasis ( phi[i], ord),
3500             if 1 < charsets_printlevel then (
3501                print("ivd: finite basis found:",length ( phi),i, qs)),
3502             for j : i + 1 thru length ( phi)  do
3503                if  not  member(j, mem) and length ( phi[i])  # length ( phi[j])
3504                    then (
3505                      if charsets_remseta ( qs, phi[j], ord)  =  setify1([])
3506                                   then (
3507                         mem :  setify1([ j, (mem) ]))),
3508             psi :  [ (psi), qs ])),
3509       (sort (  [ (psi) ], charsets_lenord)))
3510    else ( [])
3511 )$
3512 /* prime basis of cs */
3513
3514 charsets_primebasis ( cs,ord):=block(
3515    [ _is],
3516       if length ( cs)  = length ( ord)
3517           then (
3518             _is :  setify1([]))
3519
3520          else (
3521             _is : charsets_initialset ( cs, ord)),
3522       charsets_saturbasis ( cs, _is, ord)
3523 )$
3524 /* basis of saturation of ps wrt js */
3525
3526 charsets_saturbasis ( ps,js,ord):=block(
3527    [ qs,gb,zz,j],
3528    if js #  setify1([])  then (
3529       qs :  [ (ps), makelist(charsets_z[j] * js[j] - 1,j,1,length ( js) ) ],
3530 /* is this just a partial reverse */
3531       zz :  makelist(charsets_z[length ( js)  - j + 1],j,1,length ( js)),
3532       gb : poly_reduced_grobner ( qs,  [ (zz), reverse(ord) ]),
3533       qs :  [],
3534       for j : 1 thru length ( gb)  do
3535          if charsets_setdifference (  setify1([ (zz) ]), listofvars ( gb[j]) )
3536                 =  setify1([ (zz) ])
3537              then (
3538                qs :  [ gb[j], qs ]))
3539          else ( qs : ps),
3540       qs
3541 )$
3542 /*  primary decomposition of polynomial ideal generated by ps  */
3543 /*  November 1995 */
3544
3545 charsets_pridealdec ( qs,ord,medset):=block(
3546    [ ps,phi,psi,fs,pi,ph,e,i,ss,S,j,gb,gb1,urd,vrd,f,gs,k,n,_ind,indb,con,sep,wtd,tc,tco,last],
3547    ps : charsets_setdifference (  setify1([ (charsets_expand ( qs)) ]),  setify1([ 0 ]) ),
3548    ps : poly_reduced_grobner ( ps, charsets_reverse ( ord), lex),
3549    phi :  [  [ ps, 1, 0,  setify1([ 1 ])  ]  ],
3550    psi :  setify1([]),
3551    con : false,
3552    for n : 1 while phi #  [] do (
3553       fs : phi[1],
3554       phi :  [ part(phi, makelist(i,i,2,length ( phi))) ], /* rest? */
3555       sep : fs[2],
3556       wtd : fs[3],
3557       tco : fs[4],
3558       tc : tco,
3559       fs : fs[1],
3560       _ind : true,
3561       if 1 < n then (
3562          fs :  poly_reduced_grobner  ( fs, charsets_reverse ( ord), lex),
3563          if member(1, fs) then ( _ind : false)
3564          else if charsets_contain ( fs, ph, ord, plex)  then (
3565             _ind : false)),
3566       if _ind = true then (
3567          pi :  [ charsets_irrvardec ( fs, ord, medset, 0, 0)  ],
3568          if pi #  [  []  ]  then (
3569             e : length ( pi),
3570             if 1 < charsets_printlevel then (
3571                print("pid: ivd finished at",time (),e)),
3572             for i : 1 thru e do (
3573                ss :  setify1([]),
3574                S :  setify1([]),
3575                if 1 < e then (
3576                   for j : 1 thru e do
3577                      if j # i then (
3578                         ss : charsets_union ( ss,  setify1([ charsets_takele ( charsets_setdifference (  setify1([ (pi[j]) ]),  setify1([ (pi[i]) ]) ), pi[i], ord)  ]) )),
3579                   ss : charsets_prod ( charsets_factorps ( ss) ),
3580                   gb : charsets_saturbasis ( fs,  setify1([ ss ]), ord))
3581                else (
3582                   ss : 1,
3583                   gb : poly_reduced_grobner ( fs, charsets_reverse ( ord), lex)),
3584                   urd : charsets_maxindset ( pi[i], ord),
3585                   if urd =  []  then ( f : 1)
3586                   else (
3587                      vrd : (charsets_setdifference (  setify1([ (ord) ]),  setify1([ (urd) ]) )),
3588                      gb1 : poly_reduced_grobner ( gb,  [ vrd, (urd) ], lex),
3589                      f : lcm ( map(part( grobner::leadmon ( gb1[j],  [ vrd ] ), 1)/* $ ( ( j = 1 .. length ( gb1) ))*/))),
3590                   f : charsets_prod ( charsets_factorps (  setify1([ f ]) ) ),
3591                   gs : charsets_saturbasis ( gb,  setify1([ f ]), ord),
3592                   k : charsets_exponent ( gb, f, gs, ord),
3593                   if 1 < charsets_printlevel then (
3594                      print("pid: primary component found:",length ( psi)  + 1,length ( phi),e,i),
3595                      print ( gs)),
3596                   if psi =  setify1([])  then (
3597                      psi :  setify1([  [ gs, pi[i] ]  ]), ph : gs)
3598                   else (
3599                      charsets_union ( psi,  setify1([  [ gs, pi[i] ]  ]) ),
3600                      if % # psi then (
3601                         psi : %,
3602                         ph : charsets_idealint ( ph, gs, ord),
3603                         if charsets_contain ( ps, ph, ord, plex)
3604                             then ( con : true, return(false)))),
3605                   tco : charsets_idealint ( tco, charsets_saturbasisR ( gs,  setify1([ sep
3606                             * ss ]), ord), ord),
3607                   if 0 < charsets_class ( f, ord)  then (
3608                      if charsets_class ( sep * ss, ord)  = 0 then (
3609                         charsets_union (  setify1([ (gb) ]),  setify1([ f  ^ k ]) ))
3610                   else (
3611                      charsets_saturbasis ( charsets_union (  setify1([ (gb) ]),  setify1([ f
3612                          ^ k ]) ),  setify1([ sep * ss ]), ord)),
3613                   indb :  not  charsets_contain ( %, charsets_idealintR ( ph, tc, ord), ord, plex),
3614                   if indb = true then (
3615                      indb : charsets_satisfy ( charsets_union (  setify1([ (gb) ]),  setify1([ f
3616                                         ^ k ]) ), sep * ss, ord),
3617                      if indb  then (
3618                         phi :  [ (phi),  [ charsets_union (  setify1([ (gb) ]),  setify1([ f
3619                                               ^ k ]) ), sep * ss, wtd + 2, tc ]  ]))),
3620             S : charsets_union ( S,  setify1([ ss ^ charsets_exponent ( fs, ss, gb, ord)  ]) )),
3621             if con then return(false),
3622             S : charsets_union (  setify1([ (fs) ]), S),
3623             if  not  member(1, S) then (
3624                if charsets_class ( sep, ord)  = 0 then  S
3625                else (
3626                   charsets_saturbasis ( S,  setify1([ sep ]), ord)),
3627                   not  charsets_contain ( %, charsets_idealintR ( ph, tco, ord), ord, plex),
3628                   if % then (
3629                      phi :  [ (phi),  [ S, sep, wtd + 2, tco ]  ]))),
3630          phi : sort ( phi, charsets_wtdorder))),
3631    if psi =  setify1([])  or psi =  setify1([  [ 1 ]  ])  then ( [])
3632    else (
3633       (sort (  [ charsets_remred (  [ (charsets_setdifference ( psi,  setify1([  [ 1 ]  ]) )) ], ord)  ], charsets_lenord)))
3634 )$
3635 /* weight ordering */
3636
3637 charsets_wtdorder ( a,b):=block(
3638       if a[3] < b[3]
3639           then (
3640             true)
3641
3642          else (
3643             false)
3644 )$
3645 /* radical ideal membership test - is zero(ps/g) empty? */
3646
3647 charsets_satisfy ( ps,g,ord):=block(
3648       if charsets_class ( g, ord)  = 0 then (
3649             true)
3650          else if member(g, ps) or member(g ^ 2, ps) or member(g ^ 3, ps) then (
3651             false)
3652          else if charsets_contain ( ps,  setify1([ g ]), ord, plex)  then (
3653             false)
3654          else if charsets_eicsTest (  [ ps, g ], ord, charsets_charsetn)  =  []  then (
3655             false)
3656          else (
3657             true)
3658 )$
3659 /* remove redundant poly sets from phi */
3660
3661 charsets_remred ( phi,ord):=block(
3662       (charsets_remrsub ( charsets_reverse ( charsets_remrsub ( phi, ord) ), ord))
3663 )$
3664 /* main subroutine for remred */
3665
3666 charsets_remrsub ( phi,ord):=block(
3667    [ ph,mem,i,j],
3668       ph :  [],
3669       mem :  setify1([]),
3670       for i : 1 thru length ( phi)  do
3671          if  not  member(i, mem)
3672              then (
3673                for j : i + 1 thru length ( phi)  do
3674                   if  not  member(j, mem)
3675                       then (
3676                         if charsets_contain ( phi[j][1], phi[i][1], ord, plex)
3677                             then (
3678                               mem :  setify1([ (mem), j ]))),
3679                ph :  [ (ph), phi[i] ]),
3680       ph
3681 )$
3682 /* test: does ideal(ps) contain idealqs? */
3683
3684 charsets_contain ( ps,qs,ord,pt):=block(
3685    [ X,q],
3686       if 1 < charsets_printlevel
3687           then (
3688             print("containment test starts at",time ())),
3689       X : charsets_reverse ( ord),
3690       for q in qs do
3691          if grobner::normalf ( q, ps, X, pt)  # 0
3692              then (
3693                return ( false)),
3694       true
3695 )$
3696 /* reverse a list */
3697
3698 charsets_reverse ( ord):= reverse(ord)$
3699
3700 /* take from qs an element which is not in ideal(gb) */
3701
3702 charsets_takele ( gs,gb,ord):=block(
3703    [ X,ps,i],
3704       X : charsets_reverse ( ord),
3705       ps : charsets_reverse ( sort (  [ (gs) ], charsets_nopsord) ),
3706       for i : 1 thru length ( ps)  do
3707          if grobner::normalf ( ps[i], gb, X, plex)  # 0
3708              then (
3709                return ( ps[i])),
3710       error ("Check Here!")
3711 )$
3712 /* maximal independent set modulo gb  */
3713
3714 charsets_maxindset ( gb,ord):=block(
3715    [ j,x,g,_ind,u,urd],
3716       urd :  setify1([]),
3717       for x in ord do
3718          _ind : true,
3719          u : charsets_union ( urd,  setify1([ x ]) ),
3720          for g in gb do
3721             if charsets_setdifference ( listofvars ( part( grobner::leadmon ( g, charsets_reverse ( ord), plex), 2)), u)
3722                 =  setify1([])
3723                 then (
3724                   _ind : false,
3725                   break),
3726          if _ind
3727              then (
3728                urd : u),
3729        [ (urd) ]
3730 )$
3731 /* compute the exponent */
3732
3733 charsets_exponent ( ps,f,qs,ord):=block(
3734    [ k],
3735       for k : 1 do
3736          if 50 < k
3737              then (
3738                print("Check PID exponent!",k)),
3739          if  setify1([ (charsets_idealquo ( ps, f ^ k, ord)) ])  =  setify1([ qs ])
3740              then (
3741                break),
3742       if 2 < charsets_printlevel
3743           then (
3744             print("pid: exponent found:",k)),
3745       k
3746 )$
3747 /* ideal quotient  */
3748
3749 charsets_idealquo ( ps,f,ord):=block(
3750    [ qs,j,last],
3751       if  charsets_operatorp(f, charsets_ListOrSet)
3752           then (
3753             qs : charsets_idealquo ( ps, f[1], ord),
3754             if 1 < length ( f)
3755                 then (
3756                    [ part(f, makelist(i,i,2,length ( f))) ],
3757                   qs : charsets_idealint ( qs, charsets_idealquo ( ps, %, ord), ord)),
3758             qs)
3759          else (
3760             qs : charsets_idealint ( ps, f, ord),
3761             poly_reduced_grobner (  setify1([ map(simplify ( qs[j] / f) /* $ ( ( j = 1 .. length ( qs) ))*/) ]), charsets_reverse ( ord), lex))
3762 )$
3763 /* ideal charsets_intersection */
3764
3765 charsets_idealint ( ps,f,ord):=block(
3766    [ qs,gb,zz,j],
3767       if 1 < charsets_printlevel
3768           then (
3769             print("ideal charsets_intersection starts at",time ())),
3770       if  charsets_operatorp(f, charsets_ListOrSet)
3771           then (
3772             qs :  [ map(zz * ps[j]/* $ ( ( j = 1 .. length ( ps) ))*/), map(( ( 1 - zz))
3773                 * f[j]/* $ ( ( j = 1 .. length ( f) ))*/) ])
3774
3775          else (
3776             qs :  [ map(zz * ps[j]/* $ ( ( j = 1 .. length ( ps) ))*/), ( ( 1 - zz)) * f ]),
3777       gb : poly_reduced_grobner ( qs,  [ zz, map(ord[length ( ord)  - j + 1]/* $ ( ( j = 1 .. length ( ord) ))*/) ], lex),
3778       qs :  [],
3779       for j : 1 thru length ( gb)  do
3780          if charsets_setdifference (  setify1([ zz ]), listofvars ( gb[j]) )  =  setify1([ zz ])
3781              then (
3782                qs :  [ gb[j], qs ]),
3783       if 1 < charsets_printlevel
3784           then (
3785             print("ideal charsets_intersection finished at",time ())),
3786       qs
3787 )$
3788 /* ideal charsets_intersection - with remember */
3789 /* I think this should be memoized */
3790 charsets_idealintR ( ps,f,ord):=block(
3791    [ qs,gb,zz,j],
3792       if 1 < charsets_printlevel
3793           then (
3794             print("ideal charsets_intersection starts at",time ())),
3795       if  charsets_operatorp(f, charsets_ListOrSet)
3796           then (
3797             qs :  [ map(zz * ps[j]/* $ ( ( j = 1 .. length ( ps) ))*/), map(( ( 1 - zz))
3798                 * f[j]/* $ ( ( j = 1 .. length ( f) ))*/) ])
3799
3800          else (
3801             qs :  [ map(zz * ps[j]/* $ ( ( j = 1 .. length ( ps) ))*/), ( ( 1 - zz)) * f ]),
3802       gb : poly_reduced_grobner ( qs,  [ zz, map(ord[length ( ord)  - j + 1]/* $ ( ( j = 1 .. length ( ord) ))*/) ], lex),
3803       qs :  [],
3804       for j : 1 thru length ( gb)  do
3805          if charsets_setdifference (  setify1([ zz ]), listofvars ( gb[j]) )  =  setify1([ zz ])
3806              then (
3807                qs :  [ gb[j], qs ]),
3808       if 1 < charsets_printlevel
3809           then (
3810             print("ideal charsets_intersection finished at",time ())),
3811       qs
3812 )$
3813 /* basis of the saturation of ideal(ps) wrt js */
3814 /* Another memoized function */
3815 charsets_saturbasisR ( ps,js,ord):=block(
3816    [ qs,gb,zz,j],
3817    if js #  setify1([]) then (
3818       qs :  [ (ps), map(charsets_z[j] * js[j] - 1, j,1,length ( js) )  ],
3819       zz :  makelist(charsets_z[ length ( js)  - j + 1], j, 1,length ( js) ),
3820       gb : poly_reduced_grobner ( qs,  [ (zz), map(ord[length ( ord)
3821                 - j + 1]/* $ ( ( j = 1 .. length ( ord) ))*/) ], lex),
3822       qs :  [],
3823       for j : 1 thru length ( gb)  do(
3824          if charsets_setdifference (  setify1([ (zz) ]), listofvars ( gb[j]) )
3825                 =  setify1([ (zz) ])
3826              then (
3827                qs :  [ gb[j], qs ])))
3828       else (
3829             qs : ps),
3830       qs
3831 )$
3832 /* the extended irreducible char series of polyset ps */
3833 /* test: is zero(ps) empty? */
3834
3835 charsets_eicsTest ( ps,ord,medset):=block(
3836    [ qs,cs,__is,iss,n,i,j,qhi,qsi,r,rr,factorset,mind,fset,_ind,ts,den,last],
3837    if  charsets_operatorp(ps[1], charsets_ListOrSet) then (
3838       qhi :  setify1([ ps ]))
3839    else ( qhi :  setify1([  [ ps, 1 ]  ])),
3840    if medset = charsets_basset then ( mind : true) else ( mind : false),
3841    qsi :  setify1([]),
3842    if 1 < charsets_printlevel then (
3843       print("radical membership test starts at",time ())),
3844    for n : 0 while qhi #  setify1([])  and qsi =  setify1([]) do (
3845       qs : qhi[1][1],
3846       if  not  mind then (
3847          if n < 20 then (
3848             last:apply(eval_string(concat("charsets_f", substring(medset,9,length ( medset)  - 1))),[ qs, ord,  [  setify1([]), listofvars ( qs)  ], fset ]),
3849             ml2(cs,factorset, charsets_removecont ( last, ord)),
3850             factorset : charsets_union ( factorset, fset[1]))
3851          else (
3852             last:apply(eval_string(concat("charsets_", substring(medset,9,length ( medset)  - 1))), [ qs, ord ]),
3853             ml2(cs,factorset, charsets_removecont ( last, ord))))
3854       else (
3855          if n < 20 then (
3856             cs : charsets_fcharseta (  [ qs ], ord, medset),
3857             factorset : part( cs[2], 2),
3858             ml2(cs,ts, charsets_removecont ( cs[1],ord)),
3859             factorset : charsets_union ( factorset, ts))
3860          else (
3861             last:charsets_charseta (  [ qs ], ord, medset),
3862             ml2(cs,factorset, charsets_removecont ( last, ord))),
3863          if 1 < charsets_printlevel then (
3864                print("characteristic set produced"),
3865                print ( cs))),
3866       if 0 < charsets_class ( cs[1], ord)  then (
3867          ml2(ts,den, charsets_irras ( cs, ord, _ind)),
3868          break("irras5"),
3869          if ts[2] = 0 then (
3870             if  not  mind then (
3871                if  not  charsets_subset ( cs, qs)  then (
3872                   cs : charsets_charseta (  setify1([ cs, qs ]), ord, medset)),
3873                if 1 < charsets_printlevel then (
3874                   print("characteristic set produced"),
3875                   print ( cs))),
3876                if 0 < charsets_class ( cs[1], ord)  then (
3877                   ml2(ts,den, charsets_irras ( cs, ord, _ind)),
3878          break("irras6"),
3879                   if ts[2] = 0 then (
3880                      __is : charsets_union ( charsets_initialset ( cs, ord), charsets_factorps ( factorset) ),
3881                      if length ( cs)  = length ( ord)  then (
3882                         rr : charsets_nopower ( qhi[1][2]))
3883                      else (
3884                         rr : charsets_nopower ( charsets_prod (  setify1([ qhi[1][2], (__is) ]) ) )),
3885                      last:charsets_premas ( rr, cs, ord),
3886                      r : charsets_simp ( last, cs, ord),
3887                      if r # 0 then (
3888                         if r = 1 then (
3889                            qsi :  setify1([ cs, (qsi) ]))
3890                            else (
3891                               qsi :  setify1([  [ cs, charsets_simpb ( r, rr)  ], (qsi) ])))))
3892                   else (
3893                      __is : charsets_factorps ( factorset),
3894                      ts :  [ 1, 0 ])),
3895          if ts[2] # 0 then (
3896                if 1 < ts[2] then (
3897                      __is : charsets_union ( charsets_initialset (  setify1([ part(cs, makelist(i,i,1,ts[2]
3898                          - 1)) ]), ord), charsets_factorps ( factorset) ))
3899                   else (
3900                      __is : charsets_factorps ( factorset))))
3901       else (
3902          __is : charsets_factorps ( factorset),
3903          ts :  [ 1, 0 ]),
3904    iss :  setify1([]),
3905    if length ( ord)  <= length ( ps)  + 1 then (
3906          for i in __is do
3907             iss :  setify1([ iss,  [  setify1([ i, qs ]), qhi[1][2] ]  ]))
3908       else (
3909          for i : 1 thru length ( __is)  do(
3910             if i = 1 then ( 1) else (
3911                   product ( __is[j], j, 1, i - 1)),
3912             iss :  setify1([ iss,  [  setify1([ qs, __is[i] ]), %
3913                 * qhi[1][2] ]  ])),
3914    if ts[2] # 0 and ts[1] #  setify1([])  then (
3915          if  not  mind then (
3916                if  not  charsets_subset ( cs, qs)  then (
3917                      charsets_charseta (  setify1([ cs, qs ]), ord, medset))
3918                   else (
3919                      cs),
3920                if % # cs then (
3921                      if ts[2] = 1 then ( cs : qs)
3922                         else (
3923                            cs :  setify1([ qs, part(cs, makelist(i,i,1,ts[2]
3924                                - 1)) ])))),
3925          for i in ts[1] do
3926             iss :  setify1([ iss,  [  setify1([ cs, i ]), charsets_prod (  setify1([ den, qhi[1][2], (__is) ]) )  ]  ]),
3927          if 0 < charsets_class ( den, ord)  and ts[2] < charsets_class ( ts[1][1], ord)  then (
3928                iss :  setify1([ iss,  [  setify1([ cs, den ]), qhi[1][2] ]  ]))),
3929    if 1 < length ( qhi)  then (
3930          qhi :  setify1([ iss, part(qhi, makelist(i,i,2,length ( qhi))) ]))
3931       else (
3932          qhi : iss))),
3933    if qsi #  setify1([])
3934     then (
3935       maplist(lambda([x],x),qsi))
3936    else ( [])
3937 )$
3938 /* ordering wrt length */
3939
3940 charsets_lenord ( a,b):=block(
3941   [len_a : charsets_length(a),
3942    len_b : charsets_length(b)],
3943    if len_b  < len_a then ( true) else ( false)
3944 )$
3945
3946
3947 /* if elements are lists sort on the length of the first element */
3948 charsets_carorderp(a,b):=block([a1,b1],
3949    a1:if listp(a) then first(a) else a,
3950    b1:if listp(b) then first(b) else b,
3951    charsets_setorderp(a1,b1)
3952 )$
3953
3954 charsets_complexityp(a,b):=block(
3955    [ca:second(complexity(a)),cb:second(complexity(b))],
3956    if ca < cb then true else (
3957       if ca > cb then false
3958       else ordergreatp(a,b))
3959 )$
3960
3961 charsets_setorder_reversep(a,b):=
3962    if not mapatom(a) and not mapatom(b)
3963    then
3964      if length(a) > length(b) then true
3965      else(if length(a) < length(b) then false else ordergreatp(a,b))
3966    else
3967      ordergreatp (a, b) $
3968
3969 charsets_setorderp(a,b):=
3970    if not mapatom(a) and not mapatom(b)
3971    then
3972      if length(a) < length(b) then true
3973      else(if length(a) > length(b) then false else orderlessp(a,b))
3974    else
3975      orderlessp (a, b) $
3976
3977 charsets_canonlt:charsets_setorderp$
3978 utma():=apply(untrace,trace)$
3979 tracesets():=trace(intersect, charsets_union, charsets_complement, charsets_setdifference,charsets_symmdifference,
3980 charsets_powerset, charsets_setify, charsets_subset, charsets_setp, charsets_subsetp, charsets_disjointp,charsets_canonlt)$
3981
3982 tma():=trace(
3983 charsets_mcharset, charsets_iniset, charsets_irrvardec,
3984 charsets_trisetc, charsets_simp, charsets_reordera, charsets_reorderb,
3985 charsets_fsubtriset, charsets_newfactoras,
3986 charsets_irras, charsets_fcharser, charsets_irrcharser,
3987 charsets_fcharseta, charsets_initialset,
3988 charsets_adjoina, charsets_verify, charsets_trisersub,
3989 charsets_mcs,
3990 charsets_fcharsetsub, charsets_fcnorm,charsets_trivial,
3991 charsets_fcnormal, charsets_fcharsetn,
3992 charsets_subtriset, charsets_qirrcharser,
3993 charsets_charser, charsets_charset,
3994 charsets_subtrisetc,
3995 charsets_exirrcharser, charsets_eics, charsets_gcdex, charsets_ics,
3996 charsets_charseta, charsets_remset, charsets_mecs, charsets_qbasset,
3997 charsets_qirras, charsets_index, charsets_solveas, charsets_adjoin3,
3998 charsets_cfactor,
3999 charsets_irrassub,
4000 charsets_fwcharsetn, charsets_reorder,
4001 charsets_solvel, charsets_ecs,
4002 charsets_solveasr,
4003 charsets_fsubtrisetsub, charsets_wcharsetn,
4004 charsets_contractsub, charsets_union, charsets_summ,
4005 charsets_triser, charsets_select, charsets_triset, charsets_solveps,
4006 charsets_getfactor, charsets_mrank, charsets_qics,
4007 charsets_sfactor,
4008 charsets_fqcharsetn, charsets_ivd, charsets_charseries,
4009 charsets_newfactorassub, charsets_wbasset, charsets_ftriset,
4010 charsets_coeff, charsets_measure, charsets_simpa,
4011 charsets_simpb, charsets_fexcharser,
4012 /*
4013 charsets_brank,charsets_degreel,charsets_deg0,charsets_degpsmax,
4014 charsets_degpsmin,charsets_degord,
4015 charsets_prem1,charsets_compress,
4016 charsets_lvar,charsets_nopower,
4017 charsets_prem,charsets_initialset1,
4018 charsets_class,charsets_removecont,
4019 charsets_subrank,charsets_factorps,  charsets_movefactor,
4020 charsets_rank,charsets_removefactor,
4021 charsets_sort,charsets_initial, charsets_remseta,
4022 content,primpart,charsets_nopsord, charsets_deg1,
4023 charsets_unigcd, charsets_charsetn, charsets_basset, charsets_minus,
4024 charsets_gcd, charsets_getvars, charsets_pfactor,charsets_premas,
4025 */
4026 charsets_movefactorps,
4027 charsets_algcd,
4028 charsets_csolve,
4029 charsets_premasb,
4030 charsets_isirr,charsets_qfactor,
4031 charsets_cfactorsub,
4032 charsets_remsetb,charsets_tsfactor,charsets_qqfactor,
4033 charsets_contract, charsets_factoras,
4034 charsets_ftrisetc, charsets_qcharsetn,
4035 charsets_trank, charsets_prod, charsets_linas,charsets_primebasis,
4036 charsets_excharser, charsets_getfact,charsets_expand,
4037 charsets_equal,
4038 charsets_subset,
4039 setequal,setequal_member
4040 )$
4041 traceall():=trace(
4042 charsets_testsub,
4043 charsets_equal,
4044 charsets_remset,
4045 charsets_charset,
4046 charsets_expand,
4047 charsets_mcharset,
4048 charsets_iniset,
4049 charsets_charser,
4050 charsets_mcs,
4051 charsets_ecs,
4052 charsets_mecs,
4053 charsets_ics,
4054 charsets_eics,
4055 charsets_qics,
4056 charsets_cfactor,
4057 charsets_triser,
4058 charsets_csolve,
4059 charsets_ivd,
4060 charsets_pid,
4061 charsets_class,
4062 charsets_lvar,
4063 charsets_index,
4064 charsets_terms,
4065 charsets_initial,
4066 charsets_mrank,
4067 charsets_rank,
4068 charsets_subrank,
4069 charsets_trank,
4070 charsets_prem,
4071 charsets_prem1,
4072 charsets_premas,
4073 charsets_remseta,
4074 charsets_premasb,
4075 charsets_remsetb,
4076 charsets_reorder,
4077 charsets_reordera,
4078 charsets_reorderb,
4079 charsets_degord,
4080 charsets_degpsmax,
4081 charsets_degpsmin,
4082 charsets_deg0,
4083 charsets_deg1,
4084 charsets_sort,
4085 charsets_minus,
4086 charsets_union,
4087 charsets_prod,
4088 charsets_degreel,
4089 charsets_basset,
4090 charsets_wbasset,
4091 charsets_qbasset,
4092 charsets_brank,
4093 charsets_wbrank,
4094 charsets_qbrank,
4095 charsets_charseta,
4096 charsets_unigcd,
4097 charsets_gcd,
4098 charsets_fcharseta,
4099 charsets_fcharsetsub,
4100 charsets_prema,
4101 charsets_remsetaa,
4102 charsets_premasa,
4103 charsets_checkwith,
4104 charsets_nopower,
4105 charsets_charsetn,
4106 charsets_wcharsetn,
4107 charsets_qcharsetn,
4108 charsets_fcharsetn,
4109 charsets_fwcharsetn,
4110 charsets_fqcharsetn,
4111 charsets_triset,
4112 charsets_subtriset,
4113 charsets_movefactor,
4114 charsets_removefactor,
4115 charsets_removecont,
4116 charsets_ftriset,
4117 charsets_fsubtriset,
4118 charsets_fsubtrisetsub,
4119 charsets_trisetc,
4120 charsets_ftrisetc,
4121 charsets_subtrisetc,
4122 charsets_initialset1,
4123 charsets_compress,
4124 charsets_pfactor,
4125 charsets_sfactor,
4126 charsets_initialset,
4127 charsets_factorps,
4128 charsets_hack1,
4129 charsets_hack2,
4130 charsets_charseries,
4131 charsets_fcharser,
4132 charsets_excharser,
4133 charsets_simp,
4134 charsets_simpa,
4135 charsets_simpb,
4136 charsets_measure,
4137 charsets_fexcharser,
4138 charsets_irrcharser,
4139 charsets_subset,
4140 charsets_adjoinb,
4141 charsets_irras,
4142 charsets_irrassub,
4143 charsets_dfactors,
4144 charsets_fcnorm,
4145 charsets_fcnormal,
4146 charsets_gcdex,
4147 charsets_exirrcharser,
4148 charsets_select,
4149 charsets_adjoin3,
4150 charsets_adjoina,
4151 charsets_nopsord,
4152 charsets_contract,
4153 charsets_contractsub,
4154 charsets_has,
4155 charsets_linas,
4156 charsets_qirrcharser,
4157 charsets_qirras,
4158 charsets_cfactorsub,
4159 charsets_newfactoras,
4160 charsets_trivial,
4161 charsets_premb,
4162 charsets_newfactorassub,
4163 charsets_algcd,
4164 charsets_fcharsetna,
4165 charsets_malgcd,
4166 charsets_divide,
4167 charsets_linearas,
4168 charsets_arrange,
4169 charsets_arrangesub,
4170 charsets_movefactorps,
4171 charsets_vanish,
4172 charsets_qfactor,
4173 charsets_qqfactor,
4174 charsets_tefactor,
4175 charsets_noterms,
4176 charsets_isirr,
4177 charsets_getall,
4178 charsets_getclose,
4179 charsets_close,
4180 charsets_dividea,
4181 charsets_factoras,
4182 charsets_tsfactor,
4183 charsets_summ,
4184 charsets_coeff,
4185 charsets_getvars,
4186 charsets_solveps,
4187 charsets_verify,
4188 charsets_trisersub,
4189 charsets_solveas,
4190 charsets_solveasr,
4191 charsets_solvel,
4192 charsets_getfactor,
4193 charsets_getfact,
4194 charsets_irrvardec,
4195 charsets_primebasis,
4196 charsets_saturbasis,
4197 charsets_pridealdec,
4198 charsets_wtdorder,
4199 charsets_satisfy,
4200 charsets_remred,
4201 charsets_remrsub,
4202 charsets_contain,
4203 charsets_reverse,
4204 charsets_takele,
4205 charsets_maxindset,
4206 charsets_exponent,
4207 charsets_idealquo,
4208 charsets_idealint,
4209 charsets_idealintr,
4210 charsets_saturbasisr,
4211 charsets_eicstest,
4212 charsets_lenord,
4213 charsets_carorderp,
4214 charsets_complexityp,
4215 charsets_setorder_reversep,
4216 charsets_setorderp)$