Index: factory/libfac/charset/alg_factor.cc
===================================================================
--- factory/libfac/charset/alg_factor.cc	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/alg_factor.cc	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,794 @@
+////////////////////////////////////////////////////////////
+// emacs edit mode for this file is -*- C++ -*-
+////////////////////////////////////////////////////////////
+
+// FACTORY - Includes
+#include <factory.h>
+// Factor - Includes
+#include <tmpl_inst.h>
+#include <Factor.h>
+#include <SqrFree.h>
+#include <helpstuff.h>
+// Charset - Includes
+#include "csutil.h"
+#include "charset.h"
+#include "reorder.h"
+#include "algfactor.h"
+// some CC's need this:
+#include "alg_factor.h"
+
+void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
+
+#ifdef ALGFACTORDEBUG
+#  define DEBUGOUTPUT
+#else
+#  undef DEBUGOUTPUT
+#endif
+
+#include <libfac/factor/debug.h>
+
+static Varlist
+Var_is_in_AS(const Varlist & uord, const CFList & Astar);
+
+int getAlgVar(const CanonicalForm &f, Variable &X)
+{
+  if (f.inBaseDomain()) return 0;
+  if (f.inCoeffDomain())
+  {
+    if (f.level()!=0)
+    {
+      X= f.mvar();
+      return 1;
+    }
+    return getAlgVar(f.LC(),X);
+  }
+  if (f.inPolyDomain())
+  {
+    if (getAlgVar(f.LC(),X)) return 1;
+    for( CFIterator i=f; i.hasTerms(); i++)
+    {
+      if (getAlgVar(i.coeff(),X)) return 1;
+    }
+  }
+  return 0;
+}
+
+////////////////////////////////////////////////////////////////////////
+// This implements the algorithm of Trager for factorization of
+// (multivariate) polynomials over algebraic extensions and so called
+// function field extensions.
+////////////////////////////////////////////////////////////////////////
+
+// // missing class: IntGenerator:
+bool IntGenerator::hasItems() const
+{
+    return 1;
+}
+
+CanonicalForm IntGenerator::item() const
+//int IntGenerator::item() const
+{
+  //return current; //CanonicalForm( current );
+  return mapinto(CanonicalForm( current ));
+}
+
+void IntGenerator::next()
+{
+    current++;
+}
+
+// replacement for factory's broken psr
+static CanonicalForm
+mypsr ( const CanonicalForm &rr, const CanonicalForm &vv, const Variable & x )
+{
+  CanonicalForm r=rr, v=vv, l, test, lu, lv, t, retvalue;
+  int dr, dv, d,n=0;
+
+
+  dr = degree( r, x );
+  dv = degree( v, x );
+  if (dv <= dr) {l=LC(v,x); v = v -l*power(x,dv);}
+  else { l = 1; }
+  d= dr-dv+1;
+  while ( ( dv <= dr  ) && ( r != r.genZero()) ){
+    test = power(x,dr-dv)*v*LC(r,x);
+    if ( dr == 0 ) { r= CanonicalForm(0); }
+    else { r= r - LC(r,x)*power(x,dr); }
+    r= l*r -test;
+    dr= degree(r,x);
+    n+=1;
+  }
+  r= power(l, d-n)*r;
+  return r;
+}
+
+// replacement for factory's broken resultant
+static CanonicalForm
+resultante( const CanonicalForm & f, const CanonicalForm& g, const Variable & v )
+{
+  bool on_rational = isOn(SW_RATIONAL);
+  On(SW_RATIONAL);
+  CanonicalForm cd = bCommonDen( f );
+  CanonicalForm fz = f * cd;
+  cd = bCommonDen( g );
+  CanonicalForm gz = g * cd;
+  if (!on_rational)  Off(SW_RATIONAL);
+
+  return resultant(fz,gz,v);
+}
+
+// sqr-free routine for algebraic extensions
+// we need it! Ex.: f=c^2+2*a*c-1; as=[a^2+1]; f=(c+a)^2
+//static CFFList alg_sqrfree( const CanonicalForm & f )
+//{
+//  CFFList L;
+//
+//  L.append(CFFactor(f,1));
+//  return L;
+//}
+
+// Calculates a square free norm
+// Input: f(x, alpha) a square free polynomial over K(alpha),
+// alpha is defined by the minimal polynomial Palpha
+// K has more than S elements (S is defined in thesis; look getextension)
+static void
+sqrf_norm_sub( const CanonicalForm & f, const CanonicalForm & PPalpha,
+           CFGenerator & myrandom, CanonicalForm & s,  CanonicalForm & g,
+           CanonicalForm & R)
+{
+  Variable y=PPalpha.mvar(),vf=f.mvar();
+  CanonicalForm temp, Palpha=PPalpha, t;
+  int sqfreetest=0;
+  CFFList testlist;
+  CFFListIterator i;
+
+  DEBOUTLN(CERR, "sqrf_norm_sub:      f= ", f);
+  DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);
+  myrandom.reset();   s=f.mvar()-myrandom.item()*Palpha.mvar();   g=f;
+  R= CanonicalForm(0);
+  DEBOUTLN(CERR, "sqrf_norm_sub: myrandom s= ", s);
+
+  // Norm, resultante taken with respect to y
+  while ( !sqfreetest )
+  {
+    DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);
+    R = resultante(Palpha, g, y); R= R* bCommonDen(R);
+    DEBOUTLN(CERR, "sqrf_norm_sub: R= ", R);
+    // sqfree check ; R is a polynomial in K[x]
+    if ( getCharacteristic() == 0 )
+    {
+      temp= gcd(R, R.deriv(vf));
+      DEBOUTLN(CERR, "sqrf_norm_sub: temp= ", temp);
+      if (degree(temp,vf) != 0 || temp == temp.genZero() ){ sqfreetest= 0; }
+      else { sqfreetest= 1; }
+      DEBOUTLN(CERR, "sqrf_norm_sub: sqfreetest= ", sqfreetest);
+    }
+    else
+    {
+      DEBOUTMSG(CERR, "Starting SqrFreeTest(R)!");
+      // Look at SqrFreeTest!
+      // (z+a^5+w)^4 with z<w<a should not give sqfreetest=1 !
+      // for now we use this workaround with Factorize...
+      // ...but it should go away soon!!!!
+      Variable X;
+      if (getAlgVar(R,X))
+      {
+          testlist=factorize( R, X );
+      }
+      else
+        testlist= Factorize(R);
+      DEBOUTLN(CERR, "testlist= ", testlist);
+      testlist.removeFirst();
+      sqfreetest=1;
+      for ( i=testlist; i.hasItem(); i++)
+        if ( i.getItem().exp() > 1 && degree(i.getItem().factor(), R.mvar()) > 0) { sqfreetest=0; break; }
+      DEBOUTLN(CERR, "SqrFreeTest(R)= ", sqfreetest);
+    }
+    if ( ! sqfreetest )
+    {
+      myrandom.next();
+      DEBOUTLN(CERR, "sqrf_norm_sub generated new myrandom item: ", myrandom.item());
+      if ( getCharacteristic() == 0 ) t= CanonicalForm(mapinto(myrandom.item()));
+      else t= CanonicalForm(myrandom.item());
+      s= f.mvar()+t*Palpha.mvar(); // s defines backsubstitution
+      DEBOUTLN(CERR, "sqrf_norm_sub: testing s= ", s);
+      g= f(f.mvar()-t*Palpha.mvar(), f.mvar());
+      DEBOUTLN(CERR, "             gives g= ", g);
+    }
+  }
+}
+static void
+sqrf_agnorm_sub( const CanonicalForm & f, const CanonicalForm & PPalpha,
+           AlgExtGenerator & myrandom, CanonicalForm & s,  CanonicalForm & g,
+           CanonicalForm & R)
+{
+  Variable y=PPalpha.mvar(),vf=f.mvar();
+  CanonicalForm temp, Palpha=PPalpha, t;
+  int sqfreetest=0;
+  CFFList testlist;
+  CFFListIterator i;
+
+  DEBOUTLN(CERR, "sqrf_norm_sub:      f= ", f);
+  DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);
+  myrandom.reset();   s=f.mvar()-myrandom.item()*Palpha.mvar();   g=f;
+  R= CanonicalForm(0);
+  DEBOUTLN(CERR, "sqrf_norm_sub: myrandom s= ", s);
+
+  // Norm, resultante taken with respect to y
+  while ( !sqfreetest )
+  {
+    DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);
+    R = resultante(Palpha, g, y); R= R* bCommonDen(R);
+    DEBOUTLN(CERR, "sqrf_norm_sub: R= ", R);
+    // sqfree check ; R is a polynomial in K[x]
+    if ( getCharacteristic() == 0 )
+    {
+      temp= gcd(R, R.deriv(vf));
+      DEBOUTLN(CERR, "sqrf_norm_sub: temp= ", temp);
+      if (degree(temp,vf) != 0 || temp == temp.genZero() ){ sqfreetest= 0; }
+      else { sqfreetest= 1; }
+      DEBOUTLN(CERR, "sqrf_norm_sub: sqfreetest= ", sqfreetest);
+    }
+    else
+    {
+      DEBOUTMSG(CERR, "Starting SqrFreeTest(R)!");
+      // Look at SqrFreeTest!
+      // (z+a^5+w)^4 with z<w<a should not give sqfreetest=1 !
+      // for now we use this workaround with Factorize...
+      // ...but it should go away soon!!!!
+      Variable X;
+      if (getAlgVar(R,X))
+      {
+          testlist=factorize( R, X );
+      }
+      else
+        testlist= Factorize(R);
+      DEBOUTLN(CERR, "testlist= ", testlist);
+      testlist.removeFirst();
+      sqfreetest=1;
+      for ( i=testlist; i.hasItem(); i++)
+        if ( i.getItem().exp() > 1 && degree(i.getItem().factor(), R.mvar()) > 0) { sqfreetest=0; break; }
+      DEBOUTLN(CERR, "SqrFreeTest(R)= ", sqfreetest);
+    }
+    if ( ! sqfreetest )
+    {
+      myrandom.next();
+      DEBOUTLN(CERR, "sqrf_norm_sub generated new myrandom item: ", myrandom.item());
+      if ( getCharacteristic() == 0 ) t= CanonicalForm(mapinto(myrandom.item()));
+      else t= CanonicalForm(myrandom.item());
+      s= f.mvar()+t*Palpha.mvar(); // s defines backsubstitution
+      DEBOUTLN(CERR, "sqrf_norm_sub: testing s= ", s);
+      g= f(f.mvar()-t*Palpha.mvar(), f.mvar());
+      DEBOUTLN(CERR, "             gives g= ", g);
+    }
+  }
+}
+
+static void
+sqrf_norm( const CanonicalForm & f, const CanonicalForm & PPalpha,
+           const Variable & Extension, CanonicalForm & s,  CanonicalForm & g,
+           CanonicalForm & R)
+{
+  DEBOUTLN(CERR, "sqrf_norm:      f= ", f);
+  DEBOUTLN(CERR, "sqrf_norm: Palpha= ", PPalpha);
+  if ( getCharacteristic() == 0 )
+  {
+    IntGenerator myrandom;
+    DEBOUTMSG(CERR, "sqrf_norm: no extension, char=0");
+    sqrf_norm_sub(f,PPalpha, myrandom, s,g,R);
+    DEBOUTLN(CERR, "sqrf_norm:      f= ", f);
+    DEBOUTLN(CERR, "sqrf_norm: Palpha= ", PPalpha);
+    DEBOUTLN(CERR, "sqrf_norm:      s= ", s);
+    DEBOUTLN(CERR, "sqrf_norm:      g= ", g);
+    DEBOUTLN(CERR, "sqrf_norm:      R= ", R);
+  }
+  else if ( degree(Extension) > 0 ) // working over Extensions
+  {
+    DEBOUTLN(CERR, "sqrf_norm: degree of extension is ", degree(Extension));
+    AlgExtGenerator myrandom(Extension);
+    sqrf_agnorm_sub(f,PPalpha, myrandom, s,g,R);
+  }
+  else
+  {
+    FFGenerator myrandom;
+    DEBOUTMSG(CERR, "sqrf_norm: degree of extension is 0");
+    sqrf_norm_sub(f,PPalpha, myrandom, s,g,R);
+  }
+}
+
+static Varlist
+Var_is_in_AS(const Varlist & uord, const CFList & Astar){
+  Varlist output;
+  CanonicalForm elem;
+  Variable x;
+
+  for ( VarlistIterator i=uord; i.hasItem(); i++)
+  {
+    x=i.getItem();
+    for ( CFListIterator j=Astar; j.hasItem(); j++ )
+    {
+      elem= j.getItem();
+      if ( degree(elem,x) > 0 ) // x actually occures in Astar
+      {
+        output.append(x);
+        break;
+      }
+    }
+  }
+  return output;
+}
+
+// Look if Minimalpolynomials in Astar define seperable Extensions
+// Must be a power of p: i.e. y^{p^e}-x
+static int
+inseperable(const CFList & Astar)
+{
+  CanonicalForm elem;
+  int Counter= 1;
+
+  if ( Astar.length() == 0 ) return 0;
+  for ( CFListIterator i=Astar; i.hasItem(); i++)
+  {
+    elem= i.getItem();
+    if ( elem.deriv() == elem.genZero() ) return Counter;
+    else Counter += 1;
+  }
+  return 0;
+}
+
+// calculate gcd of f and g in char=0
+static CanonicalForm
+gcd0( CanonicalForm f, CanonicalForm g )
+{
+  int charac= getCharacteristic();
+  int save=isOn(SW_RATIONAL);
+  setCharacteristic(0); Off(SW_RATIONAL);
+  CanonicalForm ff= mapinto(f), gg= mapinto(g);
+  CanonicalForm result= gcd(ff,gg);
+  setCharacteristic(charac);
+  if (save) On(SW_RATIONAL);
+  return mapinto(result);
+}
+
+// calculate big enough extension for finite fields
+// Idea: first calculate k, such that q^k > S (->thesis, -> getextension)
+// Second, search k with gcd(k,m_i)=1, where m_i is the degree of the i'th
+// minimal polynomial. Then the minpoly f_i remains irrd. over q^k and we
+// have enough elements to plug in.
+static int
+getextension( IntList & degreelist, int n)
+{
+  int charac= getCharacteristic();
+  setCharacteristic(0); // need it for k !
+  int k=1, m=1, length=degreelist.length();
+  IntListIterator i;
+
+  for (i=degreelist; i.hasItem(); i++) m= m*i.getItem();
+  int q=charac;
+  while (q <= ((n*m)*(n*m)/2)) { k= k+1; q= q*charac;}
+  int l=0;
+  do {
+    for (i=degreelist; i.hasItem(); i++){
+      l= l+1;
+      if ( gcd0(k,i.getItem()) == 1){
+        DEBOUTLN(CERR, "getextension: gcd == 1, l=",l);
+        if ( l==length ){ setCharacteristic(charac);  return k; }
+      }
+      else { DEBOUTMSG(CERR, "getextension: Next iteration"); break; }
+    }
+    k= k+1; l=0;
+  }
+  while ( 1 );
+}
+
+// calculate a "primitive element"
+// K must have more than S elements (->thesis, -> getextension)
+static CFList
+simpleextension(const CFList & Astar, const Variable & Extension,
+                CanonicalForm & R){
+  CFList Returnlist, Bstar=Astar;
+  CanonicalForm s, g;
+
+  DEBOUTLN(CERR, "simpleextension: Astar= ", Astar);
+  DEBOUTLN(CERR, "simpleextension:     R= ", R);
+  DEBOUTLN(CERR, "simpleextension: Extension= ", Extension);
+  if ( Astar.length() == 1 ){ R= Astar.getFirst();}
+  else{
+    R=Bstar.getFirst(); Bstar.removeFirst();
+    for ( CFListIterator i=Bstar; i.hasItem(); i++){
+      DEBOUTLN(CERR, "simpleextension: f(x)= ", i.getItem());
+      DEBOUTLN(CERR, "simpleextension: P(x)= ", R);
+      sqrf_norm(i.getItem(), R, Extension, s, g, R);
+      // spielt die Repraesentation eine Rolle?
+      // muessen wir die Nachfolger aendern, wenn s != 0 ?
+      DEBOUTLN(CERR, "simpleextension: g= ", g);
+      if ( s != 0 ) DEBOUTLN(CERR, "simpleextension: s= ", s);
+      else DEBOUTLN(CERR, "simpleextension: s= ", s);
+      DEBOUTLN(CERR, "simpleextension: R= ", R);
+      Returnlist.insert(s);
+    }
+  }
+
+  return Returnlist;
+}
+
+CanonicalForm alg_lc(const CanonicalForm &f)
+{
+  if (f.level()>0)
+  {
+    return alg_lc(f.LC());
+  }
+  //assert(f.inCoeffDomain());
+  return f;
+}
+
+// the heart of the algorithm: the one from Trager
+static CFFList
+alg_factor( const CanonicalForm & f, const CFList & Astar, const Variable & vminpoly, const Varlist & oldord, const CFList & as)
+{
+  CFFList L, Factorlist;
+  CanonicalForm R, Rstar, s, g, h;
+  CFList substlist;
+
+  DEBINCLEVEL(CERR,"alg_factor");
+  DEBOUTLN(CERR, "alg_factor: f= ", f);
+
+  //out_cf("start alg_factor:",f,"\n");
+  substlist= simpleextension(Astar, vminpoly, Rstar);
+  DEBOUTLN(CERR, "alg_factor: substlist= ", substlist);
+  DEBOUTLN(CERR, "alg_factor: minpoly Rstar= ", Rstar);
+  DEBOUTLN(CERR, "alg_factor: vminpoly= ", vminpoly);
+
+  sqrf_norm(f, Rstar, vminpoly, s, g, R );
+  //out_cf("sqrf_norm R:",R,"\n");
+  //out_cf("sqrf_norm s:",s,"\n");
+  //out_cf("sqrf_norm g:",g,"\n");
+  DEBOUTLN(CERR, "alg_factor: g= ", g);
+  DEBOUTLN(CERR, "alg_factor: s= ", s);
+  DEBOUTLN(CERR, "alg_factor: R= ", R);
+  Off(SW_RATIONAL);
+  Variable X;
+  if (getAlgVar(R,X))
+  {
+    // factorize R over alg.extension with X
+//CERR << "alg: "<< X << " mipo=" << getMipo(X,Variable('X')) <<"\n";
+    if (R.isUnivariate())
+    {
+      DEBOUTLN(CERR, "alg_factor: factorize( ", R);
+      Factorlist =  factorize( R, X );
+    }
+    else
+    {
+      #if 1
+      Variable XX;
+      CanonicalForm mipo=getMipo(X,XX);
+      CFList as(mipo);
+      DEBOUTLN(CERR, "alg_factor: newfactoras( ", R);
+      int success=1;
+      Factorlist = newfactoras(R, as , success);
+      #else
+      // factor R over k
+      DEBOUTLN(CERR, "alg_factor: Factorize( ", R);
+      Factorlist = Factorize(R);
+      #endif
+    }
+  }
+  else
+  {
+    // factor R over k
+    DEBOUTLN(CERR, "alg_factor: Factorize( ", R);
+    Factorlist = Factorize(R);
+  }
+  On(SW_RATIONAL);
+  DEBOUTLN(CERR, "alg_factor: Factorize(R)= ", Factorlist);
+  if ( !Factorlist.getFirst().factor().inCoeffDomain() )
+    Factorlist.insert(CFFactor(1,1));
+  if ( Factorlist.length() == 2 && Factorlist.getLast().exp()== 1)
+  { // irreduzibel (first entry is a constant)
+    L.append(CFFactor(f,1));
+  }
+  else
+  {
+    DEBOUTLN(CERR, "alg_factor: g= ", g);
+    CanonicalForm gnew= g(s,s.mvar());
+    DEBOUTLN(CERR, "alg_factor: gnew= ", gnew);
+    g=gnew;
+    for ( CFFListIterator i=Factorlist; i.hasItem(); i++)
+    {
+      CanonicalForm fnew=i.getItem().factor();
+      fnew= fnew(s,s.mvar());
+      DEBOUTLN(CERR, "alg_factor: fnew= ", fnew);
+      DEBOUTLN(CERR, "alg_factor: substlist= ", substlist);
+      for ( CFListIterator ii=substlist; ii.hasItem(); ii++){
+        DEBOUTLN(CERR, "alg_factor: item= ", ii.getItem());
+        fnew= fnew(ii.getItem(), ii.getItem().mvar());
+        DEBOUTLN(CERR, "alg_factor: fnew= ", fnew);
+      }
+      if (degree(i.getItem().factor()) > 0 )
+      {
+          h=alg_gcd(g,fnew,as);
+        DEBOUTLN(CERR, "  alg_factor: h= ", h);
+        DEBOUTLN(CERR, "  alg_factor: oldord= ", oldord);
+        if ( degree(h) > 0 )
+        { //otherwise it's a constant
+          g= divide(g, h,as);
+          DEBOUTLN(CERR, "alg_factor: g/h= ", g);
+          DEBOUTLN(CERR, "alg_factor: s= ", s);
+          DEBOUTLN(CERR, "alg_factor: substlist= ", substlist);
+          //out_cf("new g:",g,"\n");
+          L.append(CFFactor(h,1));
+        }
+        //else printf("degree <=1\n");
+      }
+    }
+    // we are not interested in a
+    // constant (over K_r, which can be a polynomial!)
+    if (degree(g, f.mvar())>0){ L.append(CFFactor(g,1)); }
+  }
+  CFFList LL;
+  if (getCharacteristic()>0)
+  {
+    CFFListIterator i=L;
+    CanonicalForm c_fac=1;
+    CanonicalForm c;
+    for(;i.hasItem(); i++ )
+    {
+      CanonicalForm ff=i.getItem().factor();
+      c=alg_lc(ff);
+      int e=i.getItem().exp();
+      ff/=c;
+      if (!ff.isOne()) LL.append(CFFactor(ff,e));
+      while (e>0) { c_fac*=c;e--; }
+    }
+    if (!c_fac.isOne()) LL.insert(CFFactor(c_fac,1));
+  }
+  else
+  {
+    LL=L;
+  }
+  //CFFListIterator i=LL;
+  //for(;i.hasItem(); i++ )
+  //  out_cf("end alg_f:",i.getItem().factor(),"\n");
+  //printf("end alg_factor\n");
+  DEBOUTLN(CERR, "alg_factor: L= ", LL);
+  DEBDECLEVEL(CERR,"alg_factor");
+  return LL;
+}
+
+static CFFList
+endler( const CanonicalForm & f, const CFList & AS, const Varlist & uord )
+{
+  CanonicalForm F=f, g, q,r;
+  CFFList Output;
+  CFList One, Two, asnew, as=AS;
+  CFListIterator i,ii;
+  VarlistIterator j;
+  Variable vg;
+
+  for (i=as; i.hasItem(); i++)
+  {
+    g= i.getItem();
+    if (g.deriv() == 0 )
+    {
+      DEBOUTLN(CERR, "Inseperable extension detected: ", g);
+      for (j=uord; j.hasItem(); j++)
+      {
+        if ( degree(g,j.getItem()) > 0 ) vg= j.getItem();
+      }
+      // Now we have the highest transzendental in vg;
+      DEBOUTLN(CERR, "Transzendental is ", vg);
+      CanonicalForm gg=-1*g[0];
+      divrem(gg,vg,q,r); r= gg-q*vg;   gg= gg-r;
+      //DEBOUTLN(CERR, "q= ", q); DEBOUTLN(CERR, "r= ", r);
+      DEBOUTLN(CERR, "  that is ", gg);
+      DEBOUTLN(CERR, "  maps to ", g+gg);
+      One.insert(gg); Two.insert(g+gg);
+      // Now transform all remaining polys in as:
+      int x=0;
+      for (ii=i; ii.hasItem(); ii++)
+      {
+        if ( x != 0 )
+	{
+          divrem(ii.getItem(), gg, q,r);
+//          CERR << ii.getItem() << " divided by " << gg << "\n";
+          DEBOUTLN(CERR, "q= ", q); DEBOUTLN(CERR, "r= ", r);
+          ii.append(ii.getItem()+q*g); ii.remove(1);
+          DEBOUTLN(CERR, "as= ", as);
+        }
+        x+= 1;
+      }
+      // Now transform F:
+      divrem(F, gg, q,r);
+      F= F+q*g;
+      DEBOUTLN(CERR, "new F= ", F);
+    }
+    else{ asnew.append(i.getItem());  }// just the identity
+  }
+  // factor F with minimal polys given in asnew:
+  DEBOUTLN(CERR, "Factor F=  ", F);
+  DEBOUTLN(CERR, "  with as= ", asnew);
+  int success=0;
+  CFFList factorlist= newcfactor(F,asnew, success);
+  DEBOUTLN(CERR, "  gives =  ", factorlist);
+  DEBOUTLN(CERR, "One= ", One);
+  DEBOUTLN(CERR, "Two= ", Two);
+
+  // Transform back:
+  for ( CFFListIterator k=factorlist; k.hasItem(); k++)
+  {
+    CanonicalForm factor= k.getItem().factor();
+    ii=One;
+    for (i=Two; i.hasItem(); i++)
+    {
+      DEBOUTLN(CERR, "Mapping ", i.getItem());
+      DEBOUTLN(CERR, "     to ", ii.getItem());
+      DEBOUTLN(CERR, "     in ", factor);
+      divrem(factor,i.getItem(),q,r); r=factor -q*i.getItem();
+      DEBOUTLN(CERR, "q= ", q); DEBOUTLN(CERR, "r= ", r);
+      factor= ii.getItem()*q +r; //
+      ii++;
+    }
+    Output.append(CFFactor(factor,k.getItem().exp()));
+  }
+
+  return Output;
+}
+
+
+// 1) prepares data
+// 2) for char=p we distinguish 3 cases:
+//           no transcendentals, seperable and inseperable extensions
+CFFList
+newfactoras( const CanonicalForm & f, const CFList & as, int &success)
+{
+  Variable vf=f.mvar();
+  CFListIterator i;
+  CFFListIterator jj;
+  CFList reduceresult;
+  CFFList result;
+
+  success=1;
+  DEBINCLEVEL(CERR, "newfactoras");
+  DEBOUTLN(CERR, "newfactoras called with f= ", f);
+  DEBOUTLN(CERR, "               content(f)= ", content(f));
+  DEBOUTLN(CERR, "                       as= ", as);
+  DEBOUTLN(CERR, "newfactoras: cls(vf)= ", cls(vf));
+  DEBOUTLN(CERR, "newfactoras: cls(as.getLast())= ", cls(as.getLast()));
+  DEBOUTLN(CERR, "newfactoras: degree(f,vf)= ", degree(f,vf));
+
+// F1: [Test trivial cases]
+// 1) first trivial cases:
+  if ( (cls(vf) <= cls(as.getLast())) ||  degree(f,vf)<=1 ){
+// ||( (as.length()==1) && (degree(f,vf)==3) && (degree(as.getFirst()==2)) )
+    DEBDECLEVEL(CERR,"newfactoras");
+    return CFFList(CFFactor(f,1));
+  }
+
+// 2) List of variables:
+// 2a) Setup list of those polys in AS having degree(AS[i], AS[i].mvar()) > 1
+// 2b) Setup variableordering
+  CFList Astar;
+  Variable x;
+  CanonicalForm elem;
+  Varlist ord, uord,oldord;
+  for ( int ii=1; ii< level(vf) ; ii++ ) { uord.append(Variable(ii));  }
+  oldord= uord; oldord.append(vf);
+
+  for ( i=as; i.hasItem(); i++ ){
+    elem= i.getItem();
+    x= elem.mvar();
+    if ( degree(elem,x) > 1){ // otherwise it's not an extension
+      Astar.append(elem);
+      ord.append(x);
+    }
+  }
+  uord= Difference(uord,ord);
+  DEBOUTLN(CERR, "Astar is: ", Astar);
+  DEBOUTLN(CERR, "ord is: ", ord);
+  DEBOUTLN(CERR, "uord is: ", uord);
+
+// 3) second trivial cases: we already prooved irr. of f over no extensions
+  if ( Astar.length() == 0 ){
+    DEBDECLEVEL(CERR,"newfactoras");
+    return CFFList(CFFactor(f,1));
+  }
+
+// 4) Try to obtain a partial factorization using prop2 and prop3
+//    Use with caution! We have to proof these propositions first!
+  // Not yet implemented
+
+// 5) Look if elements in uord actually occure in any of the minimal
+//    polynomials. If no element of uord occures in any of the minimal
+//   polynomials, we don't have transzendentals.
+  Varlist newuord=Var_is_in_AS(uord,Astar);
+  DEBOUTLN(CERR, "newuord is: ", newuord);
+
+  CFFList Factorlist;
+  Varlist gcdord= Union(ord,newuord); gcdord.append(f.mvar());
+  // This is for now. we need alg_sqrfree implemented!
+  CanonicalForm Fgcd;
+          Fgcd= alg_gcd(f,f.deriv(),Astar);
+  if ( Fgcd == 0 ) DEBOUTMSG(CERR, "WARNING: p'th root ?");
+  if (( degree(Fgcd, f.mvar()) > 0) && (!(f.deriv().isZero())) ){
+    DEBOUTLN(CERR, "Nontrivial GCD found of ", f);
+    CanonicalForm Ggcd= divide(f, Fgcd,Astar);
+    DEBOUTLN(CERR, "  split into ", Fgcd);
+    DEBOUTLN(CERR, "         and ", Ggcd);
+    Fgcd= pp(Fgcd); Ggcd= pp(Ggcd);
+    DEBDECLEVEL(CERR,"newfactoras");
+    return myUnion(newfactoras(Fgcd,as,success) , newfactoras(Ggcd,as,success));
+  }
+  if ( getCharacteristic() > 0 ){
+
+    // First look for extension!
+    IntList degreelist;
+    Variable vminpoly;
+    for (i=Astar; i.hasItem(); i++){degreelist.append(degree(i.getItem()));}
+    int extdeg= getextension(degreelist, degree(f));
+    DEBOUTLN(CERR, "Extension needed of degree ", extdeg);
+
+    // Now the real stuff!
+    if ( newuord.length() == 0 ){ // no transzendentals
+      DEBOUTMSG(CERR, "No transzendentals!");
+      if ( extdeg > 1 ){
+        CanonicalForm MIPO= generate_mipo( extdeg, vminpoly);
+        DEBOUTLN(CERR, "Minpoly produced ", MIPO);
+        vminpoly= rootOf(MIPO);
+      }
+      Factorlist= alg_factor(f, Astar, vminpoly, oldord, as);
+      DEBDECLEVEL(CERR,"newfactoras");
+      return Factorlist;
+    }
+    else if ( inseperable(Astar) > 0 ){ // Look if extensions are seperable
+      // a) Use Endler
+      DEBOUTMSG(CERR, "Inseperable extensions! Using Endler!");
+      CFFList templist= endler(f,Astar, newuord);
+      DEBOUTLN(CERR, "Endler gives: ", templist);
+      return templist;
+    }
+    else{ // we are on the save side: Use trager
+      DEBOUTMSG(CERR, "Only seperable extensions!");
+      if (extdeg > 1 ){
+        CanonicalForm MIPO=generate_mipo(extdeg, vminpoly );
+        vminpoly= rootOf(MIPO);
+        DEBOUTLN(CERR, "Minpoly generated: ", MIPO);
+        DEBOUTLN(CERR, "vminpoly= ", vminpoly);
+        DEBOUTLN(CERR, "degree(vminpoly)= ", degree(vminpoly));
+      }
+      Factorlist= alg_factor(f, Astar, vminpoly, oldord, as);
+      DEBDECLEVEL(CERR,"newfactoras");
+      return Factorlist;
+    }
+  }
+  else{ // char=0 apply trager directly
+    DEBOUTMSG(CERR, "Char=0! Apply Trager!");
+    Variable vminpoly;
+    Factorlist= alg_factor(f, Astar, vminpoly, oldord, as);
+      DEBDECLEVEL(CERR,"newfactoras");
+      return Factorlist;
+  }
+
+  DEBDECLEVEL(CERR,"newfactoras");
+  return CFFList(CFFactor(f,1));
+}
+
+CFFList
+newcfactor(const CanonicalForm & f, const CFList & as, int & success )
+{
+  Off(SW_RATIONAL);
+  CFFList Output, output, Factors=Factorize(f);
+  On(SW_RATIONAL);
+  Factors.removeFirst();
+
+  if ( as.length() == 0 ){ success=1; return Factors;}
+  if ( cls(f) <= cls(as.getLast()) ) { success=1; return Factors;}
+
+  success=1;
+  for ( CFFListIterator i=Factors; i.hasItem(); i++ )
+  {
+    output=newfactoras(i.getItem().factor(),as, success);
+    for ( CFFListIterator j=output; j.hasItem(); j++)
+      Output = myappend(Output,CFFactor(j.getItem().factor(),j.getItem().exp()*i.getItem().exp()));
+  }
+  return Output;
+}
Index: factory/libfac/charset/alg_factor.h
===================================================================
--- factory/libfac/charset/alg_factor.h	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/alg_factor.h	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,33 @@
+////////////////////////////////////////////////////////////
+// emacs edit mode for this file is -*- C++ -*-
+////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////
+
+#ifndef INCL_NEW_ALGFACTOR_H
+#define INCL_NEW_ALGFACTOR_H
+
+#include <factory.h>
+#include <tmpl_inst.h>  // for typedef's
+
+// missing class: IntGenerator:
+class IntGenerator : public CFGenerator
+{
+private:
+    int current;
+public:
+    IntGenerator() : current(0) {}
+    ~IntGenerator() {}
+    bool hasItems() const;
+    void reset() { current = 0; }
+    CanonicalForm item() const;
+    void next();
+    void operator++ () { next(); }
+    void operator++ ( int ) { next(); }
+};
+
+CFFList newfactoras( const CanonicalForm & f, const CFList & as, int &success);
+/*BEGINPUBLIC*/
+CFFList newcfactor(const CanonicalForm & f, const CFList & as, int & success );
+/*ENDPUBLIC*/
+
+#endif /* INCL_ALGFACTOR_H */
Index: factory/libfac/charset/algfactor.cc
===================================================================
--- factory/libfac/charset/algfactor.cc	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/algfactor.cc	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,182 @@
+////////////////////////////////////////////////////////////
+// emacs edit mode for this file is -*- C++ -*-
+////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////
+// FACTORY - Includes
+#include <factory.h>
+// Factor - Includes
+#include <tmpl_inst.h>
+#include <Factor.h>
+#include <helpstuff.h>
+// Charset - Includes
+#include "csutil.h"
+#include "charset.h"
+#include "reorder.h"
+#include "alg_factor.h"
+// some CC's need this:
+#include "algfactor.h"
+
+#ifdef ALGFACTORDEBUG
+#  define DEBUGOUTPUT
+#else
+#  undef DEBUGOUTPUT
+#endif
+
+#include <libfac/factor/debug.h>
+#include "timing.h"
+TIMING_DEFINE_PRINT(newfactoras_time)
+
+int hasVar(const CanonicalForm &f, const Variable &v);
+
+static CFFList
+myminus( const CFFList & Inputlist, const CFFactor & TheFactor){
+  CFFList Outputlist ;
+  CFFactor copy;
+
+  for ( CFFListIterator i=Inputlist ; i.hasItem() ; i++ ){
+    copy = i.getItem();
+    if ( copy.factor() != TheFactor.factor() ){
+      Outputlist.append( copy );
+    }
+  }
+  return Outputlist;
+}
+
+static CFFList
+myDifference(const CFFList & Inputlist1,const CFFList & Inputlist2){
+  CFFList Outputlist=Inputlist1;
+
+  for ( CFFListIterator i=Inputlist2 ; i.hasItem() ; i++ )
+    Outputlist = myminus(Outputlist, i.getItem() );
+
+  return Outputlist;
+}
+
+
+// return 1 if as is quasilinear; 0 if not.
+static int
+isquasilinear( const CFList & as ){
+  if (as.length() <= 1) { return 1; }
+  else {
+    CFList AS=as;
+    AS.removeFirst(); // we are not interested in the first elem
+    for ( CFListIterator i=AS; i.hasItem(); i++)
+      if ( degree(i.getItem()) > 1 ) return 0;
+  }
+return 1;
+}
+
+#ifdef CHARSETNADEBUG
+#  define DEBUGOUTPUT
+#else
+#  undef DEBUGOUTPUT
+#endif
+#include <libfac/factor/debug.h>
+static CFList
+charsetnA(const CFList & AS, const CFList &PS, PremForm & Remembern, const Variable & vf ){
+  CFList QS = PS, RS = PS, Cset;
+  int nas= AS.length() +1;
+
+  DEBOUTLN(CERR, "charsetnA: called with ps= ", PS);
+  while ( ! RS.isEmpty() ) {
+    Cset = BasicSet( QS );
+    DEBOUTLN(CERR, "charsetnA: Cset= ", Cset);
+    Cset=Union(Cset,AS);
+    DEBOUTLN(CERR, "charsetnA: Cset= ", Cset);
+    Remembern.FS1 = Union(Remembern.FS1, initalset1(Cset));
+    DEBOUTLN(CERR, "charsetnA: Remembern.FS1= ", Remembern.FS1);
+    DEBOUTLN(CERR, "charsetnA: Remembern.FS2= ", Remembern.FS2);
+    RS = CFList();
+    if ( Cset.length() == nas && degree(Cset.getLast(),vf) > 0 ) {
+      CFList D = Difference( QS, Cset );
+      DEBOUT(CERR, "charsetnA: Difference( ", QS);
+      DEBOUT(CERR, " , ", Cset);
+      DEBOUTLN(CERR, " ) = ", D);
+      for ( CFListIterator i = D; i.hasItem(); ++i ) {
+        CanonicalForm r = Prem( i.getItem(), Cset );
+        DEBOUT(CERR,"charsetnA: Prem(", i.getItem()  );
+        DEBOUT(CERR, ",", Cset);
+        DEBOUTLN(CERR,") = ", r);
+        if ( r != 0 ){
+          //removefactor( r, Remembern );
+            RS=Union(RS,CFList(r));
+        }
+      }
+      if ( ! checkok(RS,Remembern.FS2)) return CFList(CanonicalForm(1));
+      DEBOUTLN(CERR, "charsetnA: RS= ", RS);
+      //QS = Union( QS, RS );
+      QS=Union(AS,RS); QS.append(Cset.getLast());
+      DEBOUTLN(CERR, "charsetnA: QS= Union(QS,RS)= ", QS);
+    }
+    else{ return CFList(CanonicalForm(1)); }
+  }
+  DEBOUTLN(CERR, "charsetnA: Removed factors: ", Remembern.FS2);
+  DEBOUTLN(CERR, "charsetnA: Remembern.FS1: ", Remembern.FS1);
+
+  return Cset;
+}
+
+#ifdef ALGFACTORDEBUG
+#  define DEBUGOUTPUT
+#else
+#  undef DEBUGOUTPUT
+#endif
+#include <libfac/factor/debug.h>
+// compute the GCD of f and g over the algebraic field having
+// adjoing ascending set as
+CanonicalForm
+algcd(const CanonicalForm & F, const CanonicalForm & g, const CFList & as, const Varlist & order){
+  CanonicalForm f=F;
+  int nas= as.length()+1; // the length the new char. set should have!
+  Variable vf=f.mvar();
+
+//   for ( VarlistIterator i=order; i.hasItem() ; i++ ){
+//     vf= i.getItem();
+//     nas--;
+//     if ( nas==0 ) break;
+//   }
+//   nas= as.length()+1;
+
+  DEBOUTLN(CERR, "algcd called with f= ", f);
+  DEBOUTLN(CERR, "                  g= ", g);
+  DEBOUTLN(CERR, "                 as= ", as);
+  DEBOUTLN(CERR, "              order= ", order);
+  DEBOUTLN(CERR, "         choosen vf= ", vf);
+
+  // check trivial case:
+  if ( degree(f, order.getLast())==0 || degree(g, order.getLast())==0)
+  {
+    DEBOUTLN(CERR, "algcd Result= ", 1);
+    return CanonicalForm(1);
+  }
+
+  CFList bs; bs.append(f); bs.append(g);
+  PremForm Remembern;
+  CFList cs=charsetnA(as,bs,Remembern,vf);
+  DEBOUTLN(CERR, "CharSetA((as,bs))= ", cs);
+
+//    for ( VarlistIterator i=order; i.hasItem() ; i++ ){
+//      DEBOUTLN(CERR, "vf= ", i.getItem());
+//      DEBOUTLN(CERR, "CharSetA((as,bs))= " , charsetnA(as,bs,Remembern,i.getItem()));
+//    }
+
+  CanonicalForm result;
+  if (cs.length()==nas)
+  {
+    result= cs.getLast();
+    CanonicalForm c=vcontent(result,Variable(1));
+    //CanonicalForm c=result.Lc();
+    result/=c;
+    for(CFListIterator i=as;i.hasItem(); i++ )
+    {
+      if(hasVar(result,i.getItem().mvar()))
+      {
+        c=vcontent(result,Variable(i.getItem().level()+1));
+        result/=c;
+      }
+    }
+  }
+  else result= CanonicalForm(1);
+  DEBOUTLN(CERR, "algcd Result= ", result);
+  return result;
+}
Index: factory/libfac/charset/algfactor.h
===================================================================
--- factory/libfac/charset/algfactor.h	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/algfactor.h	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,17 @@
+////////////////////////////////////////////////////////////
+// emacs edit mode for this file is -*- C++ -*-
+////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////
+
+#ifndef INCL_ALGFACTOR_H
+#define INCL_ALGFACTOR_H
+
+#include <factory.h>
+#include <tmpl_inst.h>  // for typedef's
+
+CanonicalForm algcd(const CanonicalForm & f, const CanonicalForm & g, const CFList & as, const Varlist & order);
+/*BEGINPUBLIC*/
+int hasVar(const CanonicalForm &f, const Variable &v);
+/*ENDPUBLIC*/
+
+#endif /* INCL_ALGFACTOR_H */
Index: factory/libfac/charset/charset.cc
===================================================================
--- factory/libfac/charset/charset.cc	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/charset.cc	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,553 @@
+////////////////////////////////////////////////////////////
+// emacs edit mode for this file is -*- C++ -*-
+/////////////////////////////////////////////////////////////
+// FACTORY - Includes
+#include <factory.h>
+#ifdef HAVE_IOSTREAM
+#define CIN std::cin
+#else
+#define CIN cin
+#endif
+// Factor - Includes
+#include <SqrFree.h>
+#include <Factor.h>
+#include <interrupt.h>
+// Charset - Includes
+#include "csutil.h"
+#include "algfactor.h"
+#include "alg_factor.h"
+// Some CC's need this:
+#include "charset.h"
+
+#ifdef BASICSETDEBUG
+#  define DEBUGOUTPUT
+#else
+#  undef DEBUGOUTPUT
+#endif
+
+#include <libfac/factor/debug.h>
+#include "timing.h"
+TIMING_DEFINE_PRINT(subfactorize_time)
+
+
+// forward declarations:
+static CFList     irras(CFList & AS, int &ja, CanonicalForm & reducible);
+
+#ifdef BASICSETDEBUG
+#  define DEBUGOUTPUT
+#else
+#  undef DEBUGOUTPUT
+#endif
+#include <libfac/factor/debug.h>
+
+// the next computes a characteristic set (a basic set in Wang's sense)
+CFList
+BasicSet( const CFList &PS )
+{
+    CFList QS = PS, BS, RS;
+    CanonicalForm b;
+    int cb;
+
+    DEBOUTLN(CERR, "BasicSet: called with ps= ", PS);
+    if ( PS.length() < 2 ) return PS;
+    while ( ! QS.isEmpty() ) {
+        b = lowestRank( QS );
+        cb = rank( b );
+        DEBOUTLN(CERR, "BasicSet: choose b  = ", b);
+        DEBOUTLN(CERR, "BasicSet: it's rank = ", cb);
+        BS=Union(CFList(b),BS);//BS.append( b );
+        if ( rank( b ) == 0 )
+            return Union(PS, CFList(b)) ; // b should be the first elem!
+        else {
+            RS = CFList();
+            // QS:= {q \in QS -{B} | q is reduced wrt b}
+            // We can process whole QS, because b is not reduced wrt. b
+            for ( CFListIterator i = QS; i.hasItem(); ++i )
+                if ( degree( i.getItem(), cb ) < degree( b ) )
+                    //RS.append( i.getItem() );
+                    RS = Union(CFList(i.getItem()),RS);
+            QS = RS;
+        }
+    }
+    DEBOUTLN(CERR, "BasicSet: returning bs= ", BS);
+    return BS;
+}
+
+int
+checkok( const CFList & PS, CFList & FS2){
+  CanonicalForm elem;
+
+  for ( CFListIterator i=PS; i.hasItem(); i++){
+    elem= i.getItem();
+    for (CFListIterator j=FS2; j.hasItem(); j++){
+      if (elem == j.getItem()){
+         //        FS2= Difference(FS2,CFList(elem));
+         return 0;
+      }
+    }
+  }
+  return 1;
+}
+
+#ifdef MCHARSETNDEBUG
+#  define DEBUGOUTPUT
+#else
+#  undef DEBUGOUTPUT
+#endif
+#include <libfac/factor/debug.h>
+
+// The modified CharSet (an extended characteristic set with certain factors
+// canceled; this is a characteristic set in Wang's sense)
+CFList
+MCharSetN( const CFList &PS, PremForm & Remembern ){
+  CFList QS = PS, RS = PS, CSet, OLDCS;
+
+  DEBOUTLN(CERR, "MCharSetN: called with ps= ", PS);
+  while ( ! RS.isEmpty() ) {
+    CSet = BasicSet( QS );
+    OLDCS=CSet;
+    DEBOUTLN(CERR, "MCharSetN: CS= ", CSet);
+//     if ( getNumVars(CSet.getFirst()) > 1 ){
+//       //CSet = removecontent(CSet, Remembern);
+// #ifdef MCHARSETNDEBUG
+//       CERR << "MCharSetN: CSet= " << CSet << "\n";
+// #endif
+//     }
+    Remembern.FS1 = Union(Remembern.FS1, initalset1(CSet));
+    DEBOUTLN(CERR, "MCharSetN: Remembern.FS1= ", Remembern.FS1);
+    DEBOUTLN(CERR, "MCharSetN: Remembern.FS2= ", Remembern.FS2);
+    RS = CFList();
+    if ( rank( CSet.getFirst() ) != 0 ) {
+      CFList D = Difference( QS, CSet );
+      DEBOUT(CERR, "MCharSetN: Difference( ", QS);
+      DEBOUT(CERR, " , ", CSet);
+      DEBOUTLN(CERR, " ) = ", D);
+//CERR << "MCharSetN: Difference( " << QS << " , " << CSet << " ) = " << D << "\n";
+      //PremForm Oldremember=Remembern;
+      //PremForm Newremember=Remembern;
+      for ( CFListIterator i = D; i.hasItem(); ++i ) {
+        CanonicalForm r = Prem( i.getItem(), CSet );
+        DEBOUT(CERR,"MCharSetN: Prem(", i.getItem()  );
+        DEBOUT(CERR, ",", CSet);
+        DEBOUTLN(CERR,") = ", r);
+//CERR << "MCharSetN: Prem("<< i.getItem() << "," << CSet << ") = " << r << "\n";
+        if ( r != 0 ){
+          //removefactor( r, Newremember );
+          removefactor( r, Remembern );
+          //Remembern.FS2 = Union(Remembern.FS2, Newremember.FS2);
+          //Newremember = Oldremember;
+          //if ( cls(r) > 0 )
+            //RS=Union(CFList(r),RS);//RS.append( r );
+            RS=Union(RS,CFList(r));
+        }
+      }
+      if ( ! checkok(RS,Remembern.FS2)) return CFList(CanonicalForm(1));
+      DEBOUTLN(CERR, "MCharSetN: RS= ", RS);
+      //QS = Union( QS, RS );
+      QS = Union(OLDCS,RS);
+      DEBOUTLN(CERR, "MCharSetN: QS= Union(QS,RS)= ", QS);
+    }
+    else{ return CFList(CanonicalForm(1)); }
+  }
+  DEBOUTLN(CERR, "MCharSetN: Removed factors: ", Remembern.FS2);
+  DEBOUTLN(CERR, "MCharSetN: Remembern.FS1: ", Remembern.FS1);
+
+  return CSet;
+}
+
+
+CFList
+mcharset( const CFList &PS, PremForm & Remembern ){
+  CFList cs= MCharSetN(PS, Remembern );
+  CFList rs= remsetb(Difference(PS,cs),cs);
+
+  DEBOUTLN(CERR, "mcharset: cs= ", cs);
+  DEBOUTLN(CERR, "mcharset: rs= ", rs);
+  if ( rs.length() > 0 )
+    cs= mcharset(Union(PS,Union(cs,rs)), Remembern);
+  return cs;
+}
+
+// the "original" extended characteristic set
+CFList
+CharSet( const CFList &PS ){
+  CFList QS = PS, RS = PS, CSet;
+
+  while ( ! RS.isEmpty() ) {
+    CSet = BasicSet( QS );
+    DEBOUTLN(CERR, "CharSet: CSet= ", CSet);
+    RS = CFList();
+    if ( rank( CSet.getFirst() ) != 0 ) {
+      CFList D = Difference( QS, CSet );
+      for ( CFListIterator i = D; i.hasItem(); ++i ) {
+        CanonicalForm r = Prem( i.getItem(), CSet );
+        if ( r != 0 )  RS=Union(CFList(r),RS);//RS.append( r );
+      }
+      QS = Union( QS, RS );
+    }
+  }
+  return CSet;
+}
+
+static CFList
+charseta( const CFList & PS ){
+  CFList QS = PS, RS = PS, CSet;
+
+  while ( ! RS.isEmpty() ) {
+    CSet = CharSet( QS );
+    RS = CFList();
+    if ( rank( CSet.getFirst() ) != 0 ) {
+      CFList D = Difference( QS, CSet );
+      for ( CFListIterator i = D; i.hasItem(); ++i ) {
+        CanonicalForm r = Prem( i.getItem(), CSet );
+        if ( r != 0 )  RS=Union(CFList(r),RS);//RS.append( r );
+      }
+      QS = Union(CSet,Union( QS, RS ));
+    }
+    else return CFList(CanonicalForm(1));
+  }
+  return CSet;
+}
+
+static bool
+contractsub( const CFList & cs1, const CFList & cs2 ){
+  CFListIterator i;
+
+  for ( i=cs1; i.hasItem(); i++ )
+    if ( Prem(i.getItem(),cs2 ) != 0 ) return false;
+  CFList is=initalset1(cs1);
+  for ( i=is; i.hasItem(); i++ )
+    if ( Prem(i.getItem(),cs2 ) == 0 ) return false;
+  return true;
+}
+
+static ListCFList
+contract( const ListCFList & cs){
+  ListCFList mem,ts;
+  CFList iitem,jitem;
+
+  if ( cs.length() < 2 ) return cs;
+
+  for ( ListCFListIterator i=cs; i.hasItem(); i++ ){
+    iitem=i.getItem();
+    if ( ! member(iitem, mem))
+      for ( ListCFListIterator j=i; j.hasItem(); j++){
+        jitem=j.getItem();
+        if ( ! same( iitem, jitem ) )
+          if ( ! member(jitem, mem))
+            if ( contractsub(iitem, jitem) ){
+              ts.append(jitem); mem.append(jitem);
+            }
+            else
+              if ( contractsub(jitem, iitem) ){
+                ts.append(iitem);
+              }
+      }
+  }
+  return Minus(cs,ts);
+}
+
+static ListCFList
+adjoin(const CFList & is, const CFList & qs, const ListCFList & qh ){
+  ListCFList iss,qhi;
+  ListCFListIterator j;
+  CFList iscopy,itt;
+  CFListIterator i;
+  CanonicalForm elem;
+  int ind, length;
+
+  for ( i=is ; i.hasItem(); i++ ){
+    elem=i.getItem();
+    if ( cls(elem) > 0 ) iscopy=Union(CFList(elem),iscopy);
+  }
+  if ( iscopy.isEmpty() ) return iss;
+  qhi = MyDifference(qh,qs);
+  length = qhi.length();
+  for ( i =iscopy; i.hasItem(); i++){
+    itt = Union(qs,CFList(i.getItem()));
+    ind = 0;
+    if ( length > 0 )
+      for ( j=qhi; j.hasItem(); j++ )
+        if ( subset(j.getItem(),itt )) ind=1;
+    if ( ind == 0 ) iss.append(itt);
+  }
+  return iss;
+}
+
+static ListCFList
+adjoinb(const CFList & is, const CFList & qs, const ListCFList & qh ,const CFList & cs){
+  ListCFList iss,qhi;
+  ListCFListIterator j;
+  CFList iscopy,itt;
+  CFListIterator i;
+  CanonicalForm elem;
+  int ind, length;
+
+  for ( i=is ; i.hasItem(); i++ ){
+    elem=i.getItem();
+    if ( cls(elem) > 0 ) iscopy=Union(CFList(elem),iscopy);
+  }
+  if ( iscopy.isEmpty() ) return iss;
+  qhi = MyDifference(qh,qs);
+  length = qhi.length();
+  for ( i =iscopy; i.hasItem(); i++){
+    itt = Union(Union(qs,CFList(i.getItem())),cs);
+    ind = 0;
+    if ( length > 0 )
+      for ( j=qhi; j.hasItem(); j++ )
+        if ( subset(j.getItem(),itt )) ind=1;
+    if ( ind == 0 ) {iss.append(itt);}
+  }
+  return iss;
+}
+
+static ListCFList
+sort( const ListCFList & list_to_sort ){
+  ListCFList output,copy=list_to_sort;
+  CFList l,qs1,elem;
+
+  l = copy.getLast(); copy.removeLast();
+  if ( copy.length() == 0 ){ return ListCFList(l); }
+  for ( ListCFListIterator i=copy ; i.hasItem(); i++ ){
+    elem = i.getItem();
+    if ( elem.length() > l.length() ) {
+      output = MyUnion( ListCFList(l), output);
+      l= elem;
+    }
+    else{ output = MyUnion(output, ListCFList(elem) ); }
+  }
+  //output = MyUnion(ListCFList(l),sort(output));
+  output = MyUnion(ListCFList(l),output);
+  return output;
+}
+
+#ifdef IRRCHARSERIESDEBUG
+#  define DEBUGOUTPUT
+#else
+#  undef DEBUGOUTPUT
+#endif
+#include <libfac/factor/debug.h>
+
+ListCFList
+IrrCharSeries( const CFList &PS, int opt ){
+  CanonicalForm reducible,reducible2;
+  CFList qs,cs,factorset,is,ts;
+  ListCFList pi,ppi,qqi,qsi,iss,qhi= ListCFList(PS);
+  int nr_of_iteration=0,ts2,highestlevel=0;
+
+//  CERR << getCharacteristic() << "\n";
+  for ( CFListIterator Ps=PS; Ps.hasItem(); Ps++ )
+    if ( level(Ps.getItem() ) > highestlevel ) highestlevel = level(Ps.getItem()) ;
+//  for ( int xx=1; xx <= highestlevel; xx++)
+//   CERR << Variable(xx) ;
+//  CERR << "\n";
+//  for ( CFListIterator Ps=PS; Ps.hasItem(); Ps++ )
+//    CERR << Ps.getItem() << ", " ;//<< "\n";
+//  CERR <<  "\n";
+  while ( ! qhi.isEmpty() ) {
+    qhi=sort(qhi);
+    DEBOUTLN(CERR, "qhi is: ", qhi);
+    qs=qhi.getFirst();
+    DEBOUTLN(CERR, "qs  is: ", qs);
+    DEBOUTLN(CERR, "ppi is: ", ppi);
+    ListCFList ppi1,ppi2;
+    select(ppi,qs.length(),ppi1,ppi2);
+    DEBOUTLN(CERR, "ppi1 is: ", ppi1);
+    DEBOUTLN(CERR, "ppi2 is: ", ppi2);
+    qqi = MyUnion(ppi2,qqi);
+    DEBOUTLN(CERR, "qqi is: ", qqi);
+    if ( nr_of_iteration == 0 ){ nr_of_iteration += 1; ppi = ListCFList(); }
+    else{ nr_of_iteration += 1; ppi = MyUnion(ListCFList(qs),ppi1); }
+    DEBOUTLN(CERR,"ppi is: ", ppi);
+    PremForm Remembern;
+    cs = MCharSetN(qs,Remembern);
+    DEBOUTLN(CERR, "cs is: ", cs);
+    DEBOUTLN(CERR, "factorset is: ", Remembern.FS2);
+    cs = removecontent(cs,Remembern);
+    factorset=Remembern.FS2;
+    DEBOUTLN(CERR, "cs (after removecontent) is: ", cs);
+    DEBOUTLN(CERR, "factorset is: ", factorset);
+    // Hier: removecontent einfuegen!!!!
+    if ( cls(cs.getFirst()) > 0 ){
+      ts = irras(cs,ts2, reducible);
+
+      // INTERRUPTHANDLER
+      if ( interrupt_handle() ) return ListCFList() ;
+      // INTERRUPTHANDLER
+
+      DEBOUTLN(CERR, "ts is: ", ts);
+      DEBOUTLN(CERR, "ts2 is: ", ts2);
+      // next is preliminary: should be ==0
+      if ( ts2 <= 0 ){ //irreducible
+        if ( ! subset(cs,qs) ){
+          DEBOUTMSG(CERR, "cs is not a subset of qs");
+          cs = charseta(Union(qs,cs));
+          DEBOUTLN(CERR, "new cs is: ", cs);
+        }
+        if ( ! member(cs,pi) ){
+          pi = MyUnion(pi, ListCFList(cs));
+          DEBOUTMSG(CERR, "cs is not a member of pi");
+          DEBOUTLN(CERR, "pi is: ", pi);
+          if ( cls(cs.getFirst()) > 0 ){
+            ts = irras(cs,ts2,reducible);
+
+            // INTERRUPTHANDLER
+            if ( interrupt_handle() ) return ListCFList() ;
+            // INTERRUPTHANDLER
+
+            DEBOUTLN(CERR, "ts is: ", ts);
+            DEBOUTLN(CERR, "ts2 is: ", ts2);
+            // next is preliminary: should be ==0
+            if ( ts2 <= 0 ){ //irreducible
+              qsi = MyUnion(qsi,ListCFList(cs));
+              DEBOUTLN(CERR, "qsi is: ", qsi);
+              if ( cs.length() == highestlevel ){
+                DEBOUTLN(CERR, "cs.length() == nops(ord) :", cs.length());
+                is = factorps(factorset);
+              }
+              else{
+                DEBOUT(CERR,"cs.length() != nops(ord) :", cs.length());
+                DEBOUTLN(CERR, "  nops(ord)= ", highestlevel);
+                is = Union(initalset1(cs),factorps(factorset));
+              }
+              DEBOUTLN(CERR, "is is: ", is);
+              iss = adjoin(is,qs,qqi);
+              DEBOUTLN(CERR, "iss is: ", iss);
+            }
+          }
+          else{ iss = adjoin(factorps(factorset),qs,qqi); }
+        }
+        else{
+          DEBOUTMSG(CERR, "cs is a member of pi");
+          iss = adjoin(factorps(factorset),qs,qqi); }
+        DEBOUTLN(CERR, "iss is: ", iss);
+        DEBOUTLN(CERR, "   factorps(factorset)= ", factorps(factorset));
+        DEBOUTLN(CERR, "   qs= ", qs);
+        DEBOUTLN(CERR, "   qqi= ", qqi);
+      }
+      // next is preliminary: should be !=0
+      if ( ts2 > 0 ){
+        is = factorps(factorset);
+        DEBOUTLN(CERR, "is is: ", is);
+        if ( ts2 > 1 ){
+          // setup cst: need it later for adjoinb
+          CFList cst;
+          for ( CFListIterator i=cs ; i.hasItem(); i++){
+            if ( i.getItem() == reducible ) { break; }
+            else { cst.append(i.getItem()); }
+          }
+          is = Union(initalset1(cst), is);
+          iss = MyUnion(adjoin(is,qs,qqi), adjoinb(ts,qs,qqi,cst));
+        }
+        else{ iss = adjoin(Union(is,ts),qs,qqi); }
+        DEBOUTLN(CERR, "iss is: ", iss);
+      }
+    }
+    else{
+      iss = adjoin(factorps(factorset),qs,qqi);
+      DEBOUTMSG(CERR, "case: cs is a constant.");
+      DEBOUTLN(CERR, "  qs = ", qs);
+      DEBOUTLN(CERR, "  qqi = ", qqi);
+      DEBOUTLN(CERR, "  iss = adjoin(factorps(factorset),qs,qqi) = ",iss);
+    }
+    if ( qhi.length() > 1 ){ qhi.removeFirst(); qhi = MyUnion(qhi,iss); }
+    else{ qhi = iss; }
+    DEBOUTLN(CERR, "iss is: ", iss);
+  }
+  if ( ! qsi.isEmpty() ){
+    DEBOUTLN(CERR, "qsi before contract= ", qsi);
+    if ( opt == 0 ){
+       return contract( qsi );
+    }
+    else { return qsi; }
+  }
+  else{ return ListCFList() ; }
+}
+
+// tests for characteristic sets
+//////////////////////////////////
+#ifdef IRRASDEBUG
+#  define DEBUGOUTPUT
+#else
+#  undef DEBUGOUTPUT
+#endif
+
+#include <libfac/factor/debug.h>
+
+static CFList
+irras( CFList & AS, int & ja, CanonicalForm & reducible)
+{
+  CFFList qs;
+  CFList ts,as;
+  CanonicalForm elem;
+  int ind=1,nr=0, success=-1;
+  CFListIterator i;
+
+  ja = 0;
+  DEBOUTLN(CERR, "irras: called with: AS= ", AS);
+  for ( i=AS; i.hasItem(); i++ )
+  {
+    elem = i.getItem();
+    nr += 1;
+    DEBOUT(CERR, "irras: factoring: ", elem);
+    if ( degree(elem) > 1 ) // linear poly's are irreduzible
+    {
+      qs = factorize(elem);
+      // remove a constant
+      if (qs.getFirst().factor().degree()==0) qs.removeFirst();
+    }
+    else
+    {
+      qs=(CFFactor(elem,1));
+    }
+    DEBOUTLN(CERR, "  = ", qs);
+    // INTERRUPTHANDLER
+    if ( interrupt_handle() ) return CFList() ;
+    // INTERRUPTHANDLER
+
+    if ( (qs.length() >= 2 ) || (qs.getFirst().exp() > 1))
+    {
+      DEBOUTLN(CERR, "irras: Setting ind=0, ja= ", nr);
+      ja=nr; ind=0; reducible= elem;
+      break;
+    }
+    //    else{ as.append(elem) ; }
+  }
+  //  CERR << "ind= " << ind << "\n";
+  if ( (ind == 1) ) //&& ( as.length() > 1) )
+  {
+    if ( irreducible(AS) )
+    { // as quasilinear? => irreducible!
+      ja = 0;
+      DEBOUTLN(CERR, "as is irreducible. as= ", AS);
+    }
+    else
+    {
+      i=AS;
+      for ( nr=1; nr< AS.length(); nr++)
+      {
+        as.append(i.getItem());
+        i++;
+        if ( degree(i.getItem()) > 1 )
+        {  // search for a non linear elem
+          elem=i.getItem();
+          //if (as.length()==1)
+          //  qs = Factorize2(elem,as.getFirst());
+          //else
+            qs= newfactoras(elem,as,success);
+          if ( qs.length() > 1 || qs.getFirst().exp() > 1 )
+          { //found elem is reducible
+            reducible=elem;
+            ja=nr+1;
+            break;
+          }
+        }
+      }
+    }
+  }
+  for ( CFFListIterator k=qs; k.hasItem();k++)
+    ts.append(k.getItem().factor());
+  return ts;
+}
+
+///////////////////////////////////////////////////////////////////////////////
Index: factory/libfac/charset/charset.h
===================================================================
--- factory/libfac/charset/charset.h	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/charset.h	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,20 @@
+////////////////////////////////////////////////////////////
+// emacs edit mode for this file is -*- C++ -*-
+////////////////////////////////////////////////////////////
+
+#ifndef INCL_CHARSET_H
+#define INCL_CHARSET_H
+
+#include <factory.h>
+#include <tmpl_inst.h>  // for typedef's
+
+CFList     BasicSet( const CFList &PS );
+CFList     CharSet( const CFList &PS );
+CFList     mcharset( const CFList &PS, PremForm & Remembern );
+int        checkok( const CFList & PS, CFList & FS2);
+/*BEGINPUBLIC*/
+CFList     MCharSetN( const CFList &PS, PremForm & Remembern );
+ListCFList IrrCharSeries( const CFList &PS, int opt=0 );
+/*ENDPUBLIC*/
+
+#endif /* INCL_CHARSET_H */
Index: factory/libfac/charset/csutil.cc
===================================================================
--- factory/libfac/charset/csutil.cc	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/csutil.cc	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,983 @@
+////////////////////////////////////////////////////////////
+// emacs edit mode for this file is -*- C++ -*-
+/////////////////////////////////////////////////////////////
+// FACTORY - Includes
+#include <factory.h>
+// Factor - Includes
+#include <tmpl_inst.h>
+#include <Factor.h>
+#include <SqrFree.h>
+#include <helpstuff.h>
+#include <homogfactor.h>
+// Charset - Includes
+#include "csutil.h"
+extern void out_cf(char *s1,const CanonicalForm &f,char *s2);
+extern CanonicalForm alg_lc(const CanonicalForm &f);
+extern int hasAlgVar(const CanonicalForm &f);
+
+static bool
+lowerRank ( const CanonicalForm & f, const CanonicalForm & g, int & ind )
+{
+  int df, dg;
+  Variable vf = f.mvar(), vg = g.mvar();
+
+  if ( f.inCoeffDomain() )
+  {
+    if ( g.inCoeffDomain() ) ind= 1;
+    return true;//( vg > vf );
+  }
+  else if ( g.inCoeffDomain() ) return false;
+  else if ( vf < vg ) return true;
+  else if ( vf == vg )
+  {
+    df = degree( f ); dg = degree( g );
+    if ( df < dg ) return true;
+    else if ( df == dg ) return lowerRank( LC( f ), LC( g ) , ind);
+    else return false;
+  }
+  return false;
+}
+
+CanonicalForm
+lowestRank( const CFList & F )
+{
+  CFListIterator i = F;
+  CanonicalForm f;
+  int ind=0;
+  if ( ! i.hasItem() )        return f;
+  f = i.getItem(); ++i;
+  while ( i.hasItem() )
+  {
+    //CERR << "comparing " << f << "  and " << i.getItem()
+    // << " == " << lowerRank( i.getItem(), f, ind ) << "\n";
+    if ( lowerRank( i.getItem(), f, ind ) )
+    {
+      if ( ind )
+      {
+        CFList Itemlist= get_Terms(i.getItem());
+        CFList Flist= get_Terms(f);
+
+        // Have to further compare number of terms!
+        //CERR << "compare terms! f= " << Flist.length() << "  item= "
+        //     << Itemlist.length() <<"\n";
+        if ( Itemlist.length() < Flist.length()) f = i.getItem();
+        ind=0;
+      }
+      else
+      {
+        f = i.getItem();
+      }
+    }
+    ++i;
+  }
+  return f;
+}
+
+// old version
+// CanonicalForm
+// prem ( const CanonicalForm &f, const CanonicalForm &g )
+// {
+//   CanonicalForm ff, gg, cg;
+//   int df, dg;
+//   bool reord;
+//   Variable vf, vg, v;
+//
+//   if ( (vf = f.mvar()) < (vg = g.mvar()) ) return f;
+//   else {
+//     if ( vf == vg ) {
+//       ff = f; gg = g;
+//       reord = false;
+//       v = vg;
+//     }
+//     else {
+//       v = Variable(level(f.mvar()) + 1);
+//       ff = swapvar(f,vg,v);
+//       gg = swapvar(g,vg,v);
+//       reord=true;
+//     }
+//     cg = ini( gg, v );
+//     dg = degree( gg, v );
+//     while ( ( df = degree( ff, v ) ) >= dg )
+//       ff = cg * ff - power( v, df-dg ) * gg * LC( ff, v );
+//     if ( reord ) {
+//       return swapvar( ff, vg, v );
+//     }
+//     else
+//       return ff;
+//   }
+// }
+
+CanonicalForm
+Prem ( const CanonicalForm &f, const CanonicalForm &g )
+{
+  CanonicalForm ff, gg, l, test, lu, lv, t, retvalue;
+  int df, dg;
+  bool reord;
+  Variable vf, vg, v;
+
+  if ( (vf = f.mvar()) < (vg = g.mvar()) ) return f;
+  else
+  {
+    if ( vf == vg )
+    {
+      ff = f; gg = g;
+      reord = false;
+      v = vg;
+    }
+    else
+    {
+      v = Variable(level(f.mvar()) + 1);
+      ff = swapvar(f,vg,v);
+      gg = swapvar(g,vg,v);
+      reord=true;
+    }
+    dg = degree( gg, v );
+    df = degree( ff, v );
+    if (dg <= df) {l=LC(gg); gg = gg -LC(gg)*power(v,dg);}
+    else { l = 1; }
+    while ( ( dg <= df  ) && ( !ff.isZero()) )
+    {
+      // CERR << "Start gcd..." << "\n";
+      test = gcd(l,LC(ff));
+      //CERR << "gcd(" << l << "," << LC(ff) << ")= " << test << "\n";
+      lu = l/test; lv = LC(ff)/test;
+      t = power(v,df-dg) * gg * lv;
+      if ( df == 0 ){ ff = ff.genZero(); }
+      else { ff = ff - LC(ff)*power(v,df); }
+      ff = lu*ff - t;
+      df = degree( ff, v );
+    }
+    if ( reord )
+    {
+      retvalue= swapvar( ff, vg, v );
+    }
+    else
+    {
+      retvalue= ff;
+    }
+    return retvalue;
+  }
+}
+
+static CanonicalForm
+Sprem ( const CanonicalForm &f, const CanonicalForm &g, CanonicalForm & m, CanonicalForm & q )
+{
+  CanonicalForm ff, gg, l, test, retvalue;
+  int df, dg,n;
+  bool reord;
+  Variable vf, vg, v;
+
+  if ( (vf = f.mvar()) < (vg = g.mvar()) )
+  {
+    m=CanonicalForm(0); q=CanonicalForm(0);
+    return f;
+  }
+  else
+  {
+    if ( vf == vg )
+    {
+      ff = f; gg = g;
+      reord = false;
+      v = vg; // == x
+    }
+    else
+    {
+      v = Variable(level(f.mvar()) + 1);
+      ff = swapvar(f,vg,v); // == r
+      gg = swapvar(g,vg,v); // == v
+      reord=true;
+    }
+    dg = degree( gg, v ); // == dv
+    df = degree( ff, v ); // == dr
+    if (dg <= df) {l=LC(gg); gg = gg -LC(gg)*power(v,dg);}
+    else { l = 1; }
+    n= 0;
+    while ( ( dg <= df  ) && ( !ff.isZero()) )
+    {
+      test= power(v,df-dg) * gg * LC(ff);
+      if ( df == 0 ){ff= ff.genZero();}
+      else {ff= ff - LC(ff)*power(v,df);}
+      ff = l*ff-test;
+      df= degree(ff,v);
+      n++;
+    }
+    if ( reord )
+    {
+      retvalue= swapvar( ff, vg, v );
+    }
+    else
+    {
+      retvalue= ff;
+    }
+    m= power(l,n);
+    if ( fdivides(g,m*f-retvalue) )
+      q= (m*f-retvalue)/g;
+    else
+      q= CanonicalForm(0);
+    return retvalue;
+  }
+}
+
+CanonicalForm
+divide( const CanonicalForm & ff, const CanonicalForm & f, const CFList & as)
+{
+  CanonicalForm r,m,q;
+
+  //out_cf("divide f=",ff,"\n");
+  //out_cf("divide g=",f,"\n");
+  if (f.inCoeffDomain())
+  {
+    bool b=!isOn(SW_RATIONAL);
+    On(SW_RATIONAL);
+    q=ff/f;
+    if (b) Off(SW_RATIONAL);
+  }
+  else
+    r= Sprem(ff,f,m,q); //result in q, ignore r,m
+  //CERR << "r= " << r << "  , m= " << m << "  , q= " << q << "\n";
+  r= Prem(q,as);
+  //CERR << "r= " << r << "\n";
+  //out_cf(" ->",r,"\n");
+  return r;
+}
+
+// This function allows as to be empty; in that case, it is equivalent
+// to the previous version (treating no variables as algebraic).
+static CanonicalForm
+myfitting( const CanonicalForm &f, const CFList &as )
+{
+ CanonicalForm rem=f;
+
+ if ( !(rem.isZero()) )
+ {
+   if ( getCharacteristic() > 0 )
+     return num((rem/lc(rem)));
+   else
+   {
+     On(SW_RATIONAL);
+     CanonicalForm temp= mapinto(rem);
+//      CERR << "temp= " << temp << "\n";
+//      CERR << "lc(temp)= " << lc(temp) << "\n";
+//      CERR << "temp/lc(temp)= " << temp/lc(temp) << "\n";
+//      CERR << "num(rem/lc(rem))= " << num(rem/lc(rem)) << "\n";
+
+     // If as is of length 1, and its polynomial is level 1, then
+     // we treat the first variable as algebraic and invert the leading
+     // coefficient where this variable is part of the coefficient domain.
+
+     if (as.length() == 1 && level(as.getFirst()) == 1)
+     {
+       CanonicalForm lcoeff = temp;
+       while (level(lcoeff) > 1)
+       {
+         lcoeff = LC(lcoeff);
+       }
+       // out_cf("myfitting: lcoeff = ", lcoeff, "\n");
+
+       CanonicalForm p = as.getFirst();
+       // out_cf("myfitting: p = ", p, "\n");
+
+       CanonicalForm unused, inverse;
+
+       extgcd(lcoeff, p, inverse, unused);
+       // out_cf("myfitting: inverse = ", inverse, "\n");
+
+       // This may leave temp with non-integral coefficients,
+       // which will be cleaned up below.
+       temp = temp * inverse;
+       // out_cf("myfitting: temp = ", temp, "\n");
+     }
+
+     temp= bCommonDen(temp/lc(temp))*(temp/lc(temp));
+     Off(SW_RATIONAL);
+     rem= mapinto(temp);
+     return rem;
+   }
+ }
+ else
+   return rem;
+}
+
+CanonicalForm
+Prem( const CanonicalForm &f, const CFList &L )
+{
+  CanonicalForm rem = f;
+  CFListIterator i = L;
+  for ( i.lastItem(); i.hasItem(); i-- )
+  {
+    //CERR << "   PREM: Prem(" << rem << "," ;
+    rem = Prem( rem, i.getItem() );
+    //CERR << "   PREM: Prem(" << rem << "," << i.getItem() << ")  = " << rem << "\n";
+  }
+  return myfitting(rem, CFList());
+}
+
+CanonicalForm
+Prem( const CanonicalForm &f, const CFList &L, const CFList &as )
+{
+  CanonicalForm rem = f;
+  CFListIterator i = L;
+  for ( i.lastItem(); i.hasItem(); i-- )
+  {
+    //CERR << "   PREM: Prem(" << rem << "," ;
+    rem = Prem( rem, i.getItem() );
+    //CERR << "   PREM: Prem(" << rem << "," << i.getItem() << ")  = " << rem << "\n";
+  }
+  return myfitting(rem, as);
+}
+
+CFList
+Prem( const CFList &AS, const CFList &L )
+{
+  CFList Output;
+
+  for ( CFListIterator i=AS; i.hasItem(); i++ )
+    Output = Union(CFList(Prem(i.getItem(),L)), Output);
+
+  return Output;
+}
+
+static CanonicalForm
+premasb( const CanonicalForm & f, const CFList & as)
+{
+  CanonicalForm remd=f;
+  CFList AS=as;
+
+  if ( as.length() > 1 )
+  {
+    AS.removeFirst(); // get rid of first elem
+    CanonicalForm elem;
+    while ( ! AS.isEmpty() )
+    { // thats true for at least the first iteration
+      elem= AS.getLast();
+      remd= Prem(remd,elem);
+      AS.removeLast();
+    }
+  }
+  CanonicalForm a,b;
+  if ( mydivremt(remd, as.getFirst(), a,b )){ remd= remd.genZero();}
+  else { remd= Prem(remd, as.getFirst()); }
+
+  return remd;
+}
+
+CFList
+remsetb( const CFList & ps, const CFList & as)
+{
+  CFList output;
+  CanonicalForm elem;
+  for (CFListIterator i=ps; i.hasItem(); i++)
+  {
+    elem= premasb(i.getItem(),as);
+    if ( !elem.isZero() ) output.append(elem);
+  }
+  return output;
+}
+
+// for characteristic sets
+//////////////////////////////////
+// replace the power of factors of polys in as by 1 if any
+static CFList
+nopower( const CanonicalForm & init )
+{
+  CFFList sqrfreelist;// = Factorize(init);//SqrFree(init);
+  CFList output;
+  CanonicalForm elem;
+  int count=0;
+
+  for ( CFIterator j=init; j.hasTerms(); j++ )
+    if (!(j.coeff().isOne()) ) count++;
+  //  if ( init != 1 ){
+  //  CERR << "nopower: f is " << init << "\n";
+  //  CERR << "nopower: count is " << count << "\n";}
+  if ( count > 1 ) sqrfreelist = CFFList( CFFactor(init,1));
+  else
+  {
+    sqrfreelist = Factorize(init);
+    //sqrfreelist.removeFirst();
+  }
+  for ( CFFListIterator i=sqrfreelist; i.hasItem(); i++ )
+  {
+    elem=i.getItem().factor();
+    if ( cls(elem) > 0 ) output.append(elem);
+  }
+  return output;
+}
+
+// remove the content of polys in PS; add the removed content to
+// Remembern.FS2 ( the set of removed factors )
+CFList
+removecontent ( const CFList & PS, PremForm & Remembern )
+{
+  CFListIterator i=PS;
+  if ((!i.hasItem()) || ( cls(PS.getFirst()) == 0 )) return PS;
+
+  CFList output;
+  CanonicalForm cc,elem;
+
+  for (; i.hasItem(); i++)
+  {
+    elem = i.getItem();
+    cc = content(elem, elem.mvar());
+    if ( cls(cc) > 0 )
+    {
+      output.append(elem/cc);
+      Remembern.FS2 = Union(CFList(cc), Remembern.FS2);
+    }
+    else{ output.append(elem); }
+  }
+  return output;
+}
+
+// remove possible factors in Remember.FS1 from poly r
+// Remember.FS2 contains all factors removed before
+void
+removefactor( CanonicalForm & r , PremForm & Remembern)
+{
+  int test;
+  CanonicalForm a,b,testelem;
+  CFList testlist;
+  int n=level(r);
+  CFListIterator j ;
+
+  for ( int J=1; J<= n ; J++ )
+  {
+    testlist.append(CanonicalForm(Variable(J)));
+  }
+
+  //  testlist = Union(Remembern.FS1, testlist); // add candidates
+
+  // remove already removed factors
+  for ( j = Remembern.FS2 ; j.hasItem(); j++ )
+  {
+    testelem = j.getItem();
+    while ( 1 )
+    {
+      test = mydivremt(r,testelem,a,b);
+      if ( test && b.isZero() ) r = a;
+      else break;
+    }
+  }
+
+  // Let's look if we have other canditates to remove
+  for ( j = testlist ; j.hasItem(); j++ )
+  {
+    testelem = j.getItem();
+//    if ( testelem != r && testelem != r.mvar() ){
+    if ( testelem != r )
+    {
+      while ( 1 )
+      {
+        test = divremt(r,testelem,a,b);
+        if ( test && b.isZero() )
+        {
+          Remembern.FS2= Union(Remembern.FS2, CFList(testelem));
+          r = a;
+          if ( r == 1 ) break;
+        }
+        else break;
+      }
+    }
+  }
+  //  CERR << "Remembern.FS1 = " << Remembern.FS1 << "\n";
+  //  CERR << "Remembern.FS2 = " << Remembern.FS2 << "\n";
+  //  Remembern.FS1 = Difference(Remembern.FS1, Remembern.FS2);
+  //  CERR << "  New Remembern.FS1 = " << Remembern.FS1 << "\n";
+}
+
+
+// all irreducible nonconstant factors of a set of polynomials
+CFList
+factorps( const CFList &ps )
+{
+  CFList qs;
+  CFFList q;
+  CanonicalForm elem;
+
+  for ( CFListIterator i=ps; i. hasItem(); i++ )
+  {
+    q=Factorize(i.getItem());
+    q.removeFirst();
+    // Next can be simplified ( first (already removed) elem in q is the only constant
+    for ( CFFListIterator j=q; j.hasItem(); j++ )
+    {
+      elem = j.getItem().factor();
+      if ( getNumVars(elem) > 0 )
+        qs= Union(qs, CFList(myfitting(elem, CFList())));
+    }
+  }
+  return qs;
+}
+
+// the initial of poly f wrt to the order of the variables
+static CanonicalForm
+inital( const CanonicalForm &f )
+{
+  CanonicalForm leadcoeff;
+
+  if ( cls(f) == 0 ) {return f.genOne(); }
+  else
+  {
+    leadcoeff = LC(f,lvar(f));
+    //    if ( leadcoeff != 0 )
+    return myfitting(leadcoeff, CFList()); //num(leadcoeff/lc(leadcoeff));
+    //    else return leadcoeff;
+  }
+}
+
+// the set of all nonconstant factors of initals of polys in as
+// CFList
+// initalset(const CFList &as){
+//   CanonicalForm elem;
+//   CFList is, iss,iniset;
+
+//   for ( CFListIterator i=as ; i.hasItem(); i++ ){
+//     elem = inital(i.getItem());
+//     if ( cls(elem) > 0 ) is.append(elem);
+//   }
+//   iss = factorps(is);
+//   for ( CFListIterator j=iss; j.hasItem();j++ ){
+//     elem = j.getItem();
+//     if ( cls(elem) > 0 ) iniset.append(num(elem/lc(elem)));
+//   }
+//   return iniset;
+// }
+
+// the set of nonconstant initials of Cset
+// with certain repeated factors cancelled
+CFList
+initalset1(const CFList & Cset)
+{
+  CFList temp;
+  CFList initals;
+  CanonicalForm init;
+
+  for ( CFListIterator i = Cset ; i.hasItem(); i++ )
+  {
+    initals= nopower( inital(i.getItem()) );
+    //    init= inital(i.getItem());
+    for ( CFListIterator j = initals; j.hasItem(); j++)
+    {
+      init = j.getItem();
+      if ( cls(init) > 0 )
+        temp= Union(temp, CFList(init));
+    }
+  }
+  return temp;
+}
+
+// the set of nonconstant initials of Cset of those polys
+// not having their cls higher than reducible
+// with certain repeated factors cancelled
+CFList
+initalset2(const CFList & Cset, const CanonicalForm & reducible)
+{
+  CFList temp;
+  CFList initals;
+  CanonicalForm init;
+  int clsred = cls(reducible);
+
+  for ( CFListIterator i = Cset ; i.hasItem(); i++ )
+  {
+    init = i.getItem();
+    if ( cls(init) < clsred )
+    {
+      initals= nopower( inital(init) );
+      //    init= inital(i.getItem());
+      for ( CFListIterator j = initals; j.hasItem(); j++)
+      {
+        init = j.getItem();
+        if ( cls(init) > 0 )
+          temp= Union(temp, CFList(init));
+      }
+    }
+  }
+  return temp;
+}
+
+//replace the power of factors of poly in CF init by 1 if any
+// and return the sqrfree poly
+// static CanonicalForm
+// nopower1( const CanonicalForm & init ){
+//   CFFList returnlist=Factorize(init);
+//   CanonicalForm elem, test=init.genOne();
+//   for ( CFFListIterator i= returnlist; i.hasItem(); i++){
+//     elem = i.getItem().factor();
+//     if ( cls(elem)>0 ) test *= elem;
+//   }
+//   return test;
+// }
+
+// the sequence of distinct factors of f
+//CF pfactor( ..... )
+
+// //////////////////////////////////////////
+// // for IrrCharSeries
+
+#ifdef IRRCHARSERIESDEBUG
+#  define DEBUGOUTPUT
+#else
+#  undef DEBUGOUTPUT
+#endif
+#include <libfac/factor/debug.h>
+// examine the irreducibility of as for IrrCharSeries
+int
+irreducible( const CFList & AS)
+{
+// AS is given by AS = { A1, A2, .. Ar }, d_i = degree(Ai)
+
+// 1) we test: if d_i > 1, d_j =1 for all j<>i, then AS is irreducible.
+  bool deg1=true;
+  for ( CFListIterator i = AS ; i.hasItem(); i++ )
+  {
+    if ( degree(i.getItem()) > 1 )
+    {
+      if ( deg1 ) deg1=0;
+      else return 0; // found 2nd poly with deg > 1
+    }
+  }
+  return 1;
+}
+
+
+// select an item from PS for irras
+CFList
+select( const ListCFList & PS)
+{
+  return PS.getFirst();
+}
+
+// divide list ppi in elems having length <= and > length
+void
+select( const ListCFList & ppi, int length, ListCFList & ppi1, ListCFList & ppi2)
+{
+  CFList elem;
+  for ( ListCFListIterator i=ppi ; i.hasItem(); i++ )
+  {
+    elem = i.getItem();
+    if ( ! elem.isEmpty() )
+      if ( length <= elem.length() ){ ppi2.append(elem); }
+      else { ppi1.append(elem); }
+  }
+}
+
+
+//////////////////////////////////////////////////////////////
+// help-functions for sets
+
+// is f in F ?
+static bool
+member( const CanonicalForm &f, const CFList &F)
+{
+  for ( CFListIterator i=F; i.hasItem(); i++ )
+    if ( i.getItem() == f ) return 1;
+  return 0;
+}
+
+// are list A and B the same?
+bool
+same( const CFList &A, const CFList &B )
+{
+  CFListIterator i;
+
+  for (i = A; i.hasItem(); i++)
+    if (! member(i.getItem(), B) )  return 0;
+  for (i = B; i.hasItem(); i++)
+    if (! member(i.getItem(), A) )  return 0;
+  return 1;
+}
+
+
+// is List cs contained in List of lists pi?
+bool
+member( const CFList & cs, const ListCFList & pi )
+{
+  ListCFListIterator i;
+  CFList elem;
+
+  for ( i=pi; i.hasItem(); i++)
+  {
+    elem = i.getItem();
+    if ( same(cs,elem) ) return 1;
+  }
+  return 0;
+}
+
+// is PS a subset of Cset ?
+bool
+subset( const CFList &PS, const CFList &Cset )
+{
+  //  CERR << "subset: called with: " << PS << "   " << Cset << "\n";
+  for ( CFListIterator i=PS; i.hasItem(); i++ )
+    if ( ! member(i.getItem(), Cset) )
+    {
+      //      CERR << "subset: " << i.getItem() << "  is not a member of " << Cset << "\n";
+      return 0;
+    }
+  return 1;
+}
+
+// Union of two List of Lists
+ListCFList
+MyUnion( const ListCFList & a, const ListCFList &b )
+{
+  ListCFList output;
+  ListCFListIterator i;
+  CFList elem;
+
+  for ( i = a ; i.hasItem(); i++ )
+  {
+    elem=i.getItem();
+    if ( (! elem.isEmpty()) && ( ! member(elem,output)) )
+    {
+      output.append(elem);
+    }
+  }
+
+  for ( i = b ; i.hasItem(); i++ )
+  {
+    elem=i.getItem();
+    if ( (! elem.isEmpty()) && ( ! member(elem,output)) )
+    {
+      output.append(elem);
+    }
+  }
+  return output;
+}
+
+//if list b is member of the list of lists remove b and return the rest
+ListCFList
+MyDifference( const ListCFList & a, const CFList &b)
+{
+  ListCFList output;
+  ListCFListIterator i;
+  CFList elem;
+
+  for ( i = a ; i.hasItem(); i++ )
+  {
+    elem=i.getItem();
+    if ( (! elem.isEmpty()) && ( ! same(elem,b)) )
+    {
+      output.append(elem);
+    }
+  }
+return output;
+}
+
+// remove all elements of b from list of lists a and return the rest
+ListCFList
+Minus( const ListCFList & a, const ListCFList & b)
+{
+  ListCFList output=a;
+
+  for ( ListCFListIterator i=b; i.hasItem(); i++ )
+    output = MyDifference(output, i.getItem() );
+
+  return output;
+}
+
+#if 0
+static CanonicalForm alg_lc(const CanonicalForm &f, const CFList as)
+{
+  for(CFListIterator i=as; i.hasItem(); i++)
+  {
+    if (f.mvar()==i.getItem().mvar()) return f;
+  }
+  if (f.level()>0)
+  {
+    return alg_lc(f.LC(),as);
+  }
+  return CanonicalForm(1);
+}
+#endif
+
+CanonicalForm alg_gcd(const CanonicalForm & fff, const CanonicalForm &ggg,
+                      const CFList &as)
+{
+  CanonicalForm f=fff;
+  CanonicalForm g=ggg;
+  f=Prem(f,as);
+  g=Prem(g,as);
+  if ( f.isZero() )
+  {
+    if ( g.lc().sign() < 0 ) return -g;
+    else                     return g;
+  }
+  else  if ( g.isZero() )
+  {
+    if ( f.lc().sign() < 0 ) return -f;
+    else                     return f;
+  }
+  //out_cf("alg_gcd(",fff," , ");
+  //out_cf("",ggg,")\n");
+  CanonicalForm res;
+  // does as appear in f and g ?
+  bool has_alg_var=false;
+  for ( CFListIterator j=as;j.hasItem(); j++ )
+  {
+    Variable v=j.getItem().mvar();
+    if (hasVar(f,v)) {has_alg_var=true; /*break;*/}
+    if (hasVar(g,v)) {has_alg_var=true; /*break;*/}
+    //out_cf("as:",j.getItem(),"\n");
+  }
+  if (!has_alg_var)
+  {
+    if ((hasAlgVar(f))
+    || (hasAlgVar(g)))
+    {
+      Varlist ord;
+      for ( CFListIterator j=as;j.hasItem(); j++ )
+        ord.append(j.getItem().mvar());
+      res=algcd(f,g,as,ord);
+    }
+    else
+      res=gcd(f,g);
+    //out_cf("gcd=",res,"\n");
+    //out_cf("of f=",fff," , ");
+    //out_cf("and g=",ggg,"\n");
+
+    return res;
+  }
+
+  int mvf=f.level();
+  int mvg=g.level();
+  if (mvg > mvf)
+  {
+    CanonicalForm tmp=f; f=g; g=tmp;
+    int tmp2=mvf; mvf=mvg; mvg=tmp2;
+  }
+  if (g.inBaseDomain() || f.inBaseDomain())
+  {
+    //printf("const\n");
+    //out_cf("of f=",fff," , ");
+    //out_cf("and g=",ggg,"\n");
+    return CanonicalForm(1);
+  }
+
+  // gcd of all coefficients:
+  CFIterator i=f;
+  CanonicalForm c_gcd=i.coeff(); i++;
+  while( i.hasTerms())
+  {
+    c_gcd=alg_gcd(i.coeff(),c_gcd,as);
+    if (c_gcd.inBaseDomain()) break;
+    i++;
+  }
+  //printf("f.mvar=%d (%d), g.mvar=%d (%d)\n",f.level(),mvf,g.level(),mvg);
+  if (mvf!=mvg) // => mvf > mvg
+  {
+    res=alg_gcd(g,c_gcd,as);
+    //out_cf("alg_gcd1=",res,"\n");
+    //out_cf("of f=",fff," , ");
+    //out_cf("and g=",ggg,"\n");
+    return res;
+  }
+  // now: mvf==mvg, f.level()==g.level()
+  if (!c_gcd.inBaseDomain())
+  {
+    i=g;
+    while( i.hasTerms())
+    {
+      c_gcd=alg_gcd(i.coeff(),c_gcd,as);
+      if (c_gcd.inBaseDomain()) break;
+      i++;
+    }
+  }
+
+  //f/=c_gcd;
+  //g/=c_gcd;
+  if (!c_gcd.isOne())
+  {
+    f=divide(f,c_gcd,as);
+    g=divide(g,c_gcd,as);
+  }
+
+  CFList gg;
+  CanonicalForm r=1;
+  while (1)
+  {
+    //printf("f.mvar=%d, g.mvar=%d\n",f.level(),g.level());
+    gg=as;
+    if (!g.inCoeffDomain()) gg.append(g);
+    //out_cf("Prem(",f," , ");
+    //out_cf("",g,")\n");
+    if (g.inCoeffDomain()|| g.isZero())
+    {
+      //printf("in coeff domain:");
+      if (g.isZero()) { //printf("0\n");
+        i=f;
+        CanonicalForm f_gcd=i.coeff(); i++;
+        while( i.hasTerms())
+        {
+          f_gcd=alg_gcd(i.coeff(),f_gcd,as);
+          if (f_gcd.inBaseDomain()) break;
+          i++;
+        }
+        //out_cf("g=0 -> f:",f,"\n");
+        //out_cf("f_gcd:",f_gcd,"\n");
+        //out_cf("c_gcd:",c_gcd,"\n");
+        //f/=f_gcd;
+        f=divide(f,f_gcd,as);
+        //out_cf("f/f_gcd:",f,"\n");
+        f*=c_gcd;
+        //out_cf("f*c_gcd:",f,"\n");
+        CanonicalForm r_lc=alg_lc(f);
+        //out_cf("r_lc:",r_lc,"\n");
+        //f/=r_lc;
+        f=divide(f,r_lc,as);
+        //out_cf(" -> gcd:",f,"\n");
+        return f;
+      }
+      else { //printf("c\n");
+        return c_gcd;}
+    }
+    else if (g.level()==f.level()) r=Prem(f,gg,as);
+    else
+    {
+      //printf("diff. level: %d, %d\n", f.level(), g.level());
+      // g and f have different levels
+      //out_cf("c_gcd=",c_gcd,"\n");
+    //out_cf("of f=",fff," , ");
+    //out_cf("and g=",ggg,"\n");
+      return c_gcd;
+    }
+    //out_cf("->",r,"\n");
+    f=g;
+    g=r;
+  }
+  if (degree(f,Variable(mvg))==0)
+  {
+  // remark: mvf == mvg == f.level() ==g.level()
+    //out_cf("c_gcd=",c_gcd,"\n");
+    //out_cf("of f=",fff," , ");
+    //out_cf("and g=",ggg,"\n");
+    return c_gcd;
+  }
+  //out_cf("f=",f,"\n");
+  i=f;
+  CanonicalForm k=i.coeff(); i++;
+  //out_cf("k=",k,"\n");
+  while( i.hasTerms())
+  {
+    if (k.inCoeffDomain()) break;
+    k=alg_gcd(i.coeff(),k,as);
+    //out_cf("k=",k,"\n");
+    i++;
+  }
+  //out_cf("end-k=",k,"\n");
+  f/=k;
+  //out_cf("f/k=",f,"\n");
+  res=c_gcd*f;
+  CanonicalForm r_lc=alg_lc(res);
+  res/=r_lc;
+  //CanonicalForm r_lc=alg_lc(res,as);
+  //res/=r_lc;
+  //out_cf("alg_gcd2=",res,"\n");
+  //  out_cf("of f=",fff," , ");
+  //  out_cf("and g=",ggg,"\n");
+  //return res;
+  //if (res.level()<fff.level() && res.level() < ggg.level())
+  //  return alg_gcd(alg_gcd(fff,res,as),ggg,as);
+  //else
+    return res;
+}
Index: factory/libfac/charset/csutil.h
===================================================================
--- factory/libfac/charset/csutil.h	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/csutil.h	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,77 @@
+////////////////////////////////////////////////////////////
+// emacs edit mode for this file is -*- C++ -*-
+////////////////////////////////////////////////////////////
+
+#ifndef INCL_CSUTIL_H
+#define INCL_CSUTIL_H
+
+#include <factory.h>
+#include "algfactor.h"
+#include <tmpl_inst.h>
+
+// inline function's:
+
+inline int 
+cls ( const CanonicalForm & f ){ if ( getNumVars(f) == 0 ) {return 0;} else {return f.level();} }
+
+inline int 
+rank ( const CanonicalForm & f ){ return cls(f); }
+
+inline CanonicalForm 
+ini ( const CanonicalForm & f ){ return LC(f); }
+
+inline CanonicalForm 
+ini ( const CanonicalForm & f , const Variable v ){ return LC(f,v); }
+
+inline Variable 
+lvar ( const CanonicalForm & f ){ return f.mvar(); }
+
+///////////////////////////////////////////////////////////////
+// a class definition needed in charset.cc for Prem
+///////////////////////////////////////////////////////////////
+/*BEGINPUBLIC*/
+
+class PremForm {
+public:
+  CFList FS1;
+  CFList FS2;
+  inline PremForm& operator=( const PremForm&  value ){
+    if ( this != &value ){
+      FS1 = value.FS1;
+      FS2 = value.FS2;
+    }
+    return *this;
+  }
+};
+
+/*ENDPUBLIC*/
+
+// functions from csutil.cc:
+/*BEGINPUBLIC*/
+CanonicalForm  Prem( const CanonicalForm &f, const CanonicalForm &g );
+CanonicalForm  Prem( const CanonicalForm &f, const CFList &L );
+CFList         Prem( const CFList &AS, const CFList &L );
+CanonicalForm alg_gcd(const CanonicalForm &, const CanonicalForm &, const CFList &);
+/*ENDPUBLIC*/
+CanonicalForm  divide( const CanonicalForm & ff, const CanonicalForm & f, const CFList & as);
+CFList         remsetb( const CFList & ps, const CFList & as);
+CanonicalForm  lowestRank( const CFList & F );
+
+CFList         removecontent ( const CFList & PS, PremForm & Remembern );
+void           removefactor( CanonicalForm & r , PremForm & Remembern);
+CFList         factorps( const CFList &ps );
+//CFList         initalset(const CFList & CSet);
+CFList         initalset1(const CFList & CSet);
+CFList         initalset2(const CFList & CSet, const CanonicalForm & reducible);
+int            irreducible( const CFList & ASet);
+CFList         select( const ListCFList & PSet);
+void           select( const ListCFList & ppi, int length, ListCFList & ppi1, 
+		       ListCFList & ppi2);
+bool           same( const CFList &A, const CFList &B );
+bool           member( const CFList & cs, const ListCFList & pi );
+bool           subset( const CFList &PSet, const CFList &CSet );
+ListCFList     MyUnion( const ListCFList & a, const ListCFList &b );
+ListCFList     MyDifference( const ListCFList & a, const CFList &b);
+ListCFList     Minus( const ListCFList & a, const ListCFList &b);
+#endif /* INCL_CSUTIL_H */
+
Index: factory/libfac/charset/reorder.cc
===================================================================
--- factory/libfac/charset/reorder.cc	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/reorder.cc	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,405 @@
+////////////////////////////////////////////////////////////
+// This is really experimental. What is the best order of the variables for
+// an ideal? Clearly, if a variable only occures in one and only one polynomial
+// this variable should be assigned the highest level.
+// But what for the others?
+//
+// The strategy used here is o.k. for *most* of the examples. But there are
+// examples where this strategy doesn't give the best possible answer.
+// So, use this at your own risk! The order given may increase or decrease
+// computing time for your example!
+/////////////////////////////////////////////////////////////
+// emacs edit mode for this file is -*- C++ -*-
+////////////////////////////////////////////////////////////
+// FACTORY - Includes
+#include <factory.h>
+// Factor - Includes
+#include <tmpl_inst.h>
+#include "homogfactor.h"
+// Charset - Includes
+
+// some CC's need this:
+#include "reorder.h"
+
+#ifdef REORDERDEBUG
+#  define DEBUGOUTPUT
+#else
+#  undef DEBUGOUTPUT
+#endif
+
+#include <libfac/factor/debug.h>
+#include "timing.h"
+TIMING_DEFINE_PRINT(neworder_time)
+
+#define __ARRAY_INIT__ -1
+
+// the maximal degree of polys in PS wrt. variable x
+static int
+degpsmax( const CFList & PS, const Variable & x, 
+	  Intarray & A, Intarray & C){
+  int varlevel = level(x);
+  if ( A[varlevel] != __ARRAY_INIT__ ) return A[varlevel];
+  int max=0,temp, count=0;
+  
+  for ( CFListIterator i=PS; i.hasItem();i++){
+    temp = degree(i.getItem(),x);
+    if ( temp > max ) {max = temp; count = 0; }
+    if ( temp == max) { count += max; } // we count the number of polys
+  }
+  A[varlevel] = max; C[varlevel] = count;
+  return max;
+}
+
+// the minimal non-zero degree of polys in PS wrt. x
+// returns 0 if variable x doesn't occure in any of the polys
+static int
+degpsmin( const CFList & PS, const Variable & x, Intarray & A, Intarray & B,
+	  Intarray & C, Intarray & D){
+  int varlevel = level(x);
+  if ( B[varlevel] != __ARRAY_INIT__ ) return B[varlevel];
+  int min=degpsmax(PS,x,A,C), temp, count=0;
+
+  if (min==0) {
+    B[varlevel] = min; D[varlevel] = min;
+    return min;
+  }
+  else
+    for ( CFListIterator i=PS; i.hasItem();i++){
+      temp = degree(i.getItem(),x);
+      if ( temp < min && temp != 0 ){ min = temp; count = 0; }
+      if ( temp == min) { count += min;} // we count the number of polys
+    }
+  B[varlevel] = min; D[varlevel] = count;
+  return min;
+}
+
+// the minimal total degree of lcoeffs of polys in PS wrt. x
+// for those polys having degree degpsmin in x.
+// F will be assigned the minimal number of terms of those lcoeffs
+static int
+Tdeg( const CFList & PS, const Variable & x, Intarray & A, Intarray & B,
+      Intarray & C, Intarray & D, Intarray & E, Intarray & F){
+  int k= degpsmin(PS,x,A,B,C,D), 
+    varlevel= level(x), min=0;
+
+  if ( E[varlevel] != __ARRAY_INIT__ ) return E[varlevel];
+  if (k == 0){
+    E[varlevel]= 0; F[varlevel]= 0;
+  }
+  else{
+    int nopslc=0;
+    CFList LCdegList;
+    CanonicalForm elem;
+    CFListIterator i;
+
+    for ( i=PS; i.hasItem(); i++ ){
+      elem = i.getItem();
+      if ( degree(elem,x) == k ) LCdegList.append(LC(elem,x));
+    }
+    
+    if ( LCdegList.length() > 0 ){
+      CFList TermList;
+      int newmin, newnopslc;
+
+      min= totaldegree(LCdegList.getFirst());
+      TermList= get_Terms(LCdegList.getFirst()); 
+      nopslc= TermList.length();
+      for ( i=LCdegList; i.hasItem(); i++){
+	elem= i.getItem();
+	newmin= totaldegree(elem);
+	TermList= get_Terms(elem);
+	newnopslc= TermList.length();
+	if ( newmin < min ) min= newmin; 
+	if ( newnopslc < nopslc ) nopslc= newnopslc;
+      }
+      
+    }
+    E[varlevel]= min;
+    F[varlevel]= nopslc;
+  }
+  return min;
+}
+
+// The number of the poly in which Variable x first occures
+static int
+nr_of_poly( const CFList & PS, const Variable & x, Intarray & G){
+  int min=0, varlevel=level(x);
+  if ( G[varlevel] != __ARRAY_INIT__ ) return G[varlevel];
+  for ( CFListIterator i=PS; i.hasItem(); i++ ){
+    min += 1;
+    if ( degree(i.getItem(),x) > 0 ) break;
+  }
+  G[varlevel]= min;
+  return min;
+}
+
+static int
+degord( const Variable & x, const Variable & y, const CFList & PS,
+	Intarray & A, Intarray & B, Intarray & C, Intarray & D,
+	Intarray & E, Intarray & F, Intarray & G){
+  int xlevel= level(x), ylevel= level(y);
+
+  if      (degpsmax(PS,y,A,C) < degpsmax(PS,x,A,C))         return 1;
+  else if (degpsmax(PS,x,A,C) < degpsmax(PS,y,A,C) )        return 0;
+  else if (C[ylevel] < C[xlevel])                           return 1;
+  else if (C[xlevel] < C[ylevel])                           return 0;
+  else if (degpsmin(PS,x,A,B,C,D) < degpsmin(PS,y,A,B,C,D)) return 1;
+  else if (degpsmin(PS,y,A,B,C,D) < degpsmin(PS,x,A,B,C,D)) return 0;
+  else if (D[ylevel] < D[xlevel])                           return 1;
+  else if (D[xlevel] < D[ylevel])                           return 0;
+  else if (Tdeg(PS,y,A,B,C,D,E,F) < Tdeg(PS,x,A,B,C,D,E,F)) return 1;
+  else if (Tdeg(PS,x,A,B,C,D,E,F) < Tdeg(PS,y,A,B,C,D,E,F)) return 0;
+  else if (F[ylevel] < F[xlevel])                           return 1;
+  else if (F[xlevel] < F[ylevel])                           return 0;
+  else if (nr_of_poly(PS,x,G) <= nr_of_poly(PS,y,G))        return 1;
+  else return 0;
+
+}
+
+// determine the highest variable of all involved Variables in PS
+// NOTE:
+//    this doesn't give always the correct answer:
+//    If a variable is assigned the highest level in the definition of the
+//    original ring, but doesn't occure in any of the 
+//    polynomials, get_max_var returns the variable with a level lower than 
+//    the highest level.
+//    Is there a workaround?
+// But for the redefinition of the ring this doesn't matter due to the
+// implementation of neworder().
+
+static Variable
+get_max_var(const CFList & PS){
+  Variable x=PS.getFirst().mvar(), y;
+  for (CFListIterator i=PS; i.hasItem(); i++){
+    y = i.getItem().mvar();
+    if ( y > x ) x=y;
+  }
+  return x;
+}
+
+
+// determine if variable x is in one and only one of the polynomials in PS
+// first criterion for neworder
+static CFList
+only_in_one( const CFList & PS , const Variable & x){
+  CFList output;
+
+  for ( CFListIterator i=PS; i.hasItem(); i++ ){
+    if ( degree(i.getItem(),x) >= 1 ) output.insert(i.getItem());
+    if ( output.length() >= 2 ) break;
+  }
+  return output;
+}
+
+// initialize all Arrays (of same range) with __ARRAY_INIT__
+static void
+initArray(const int highest_level, Intarray & A, Intarray & B, Intarray & C, 
+	  Intarray & D, Intarray & E, Intarray & F, Intarray & G){ 
+
+  for ( int i=1 ; i <= highest_level; i ++){
+    A[i] = __ARRAY_INIT__; B[i] = __ARRAY_INIT__; 
+    C[i] = __ARRAY_INIT__; D[i] = __ARRAY_INIT__;
+    E[i] = __ARRAY_INIT__; F[i] = __ARRAY_INIT__;
+    G[i] = __ARRAY_INIT__;
+  }
+}
+
+// now for the second criterion; a little more complex
+//
+// idea: set up seven arrays of lenth highest_level 
+//       (of which some entries may be empty, because of only_in_one!)
+//   A saves maxdegree of Variable x in A(level(x))
+//   B saves mindegree of Variable x in B(level(x))
+//   C saves the number of polys in PS which have degree A(level(x)) in 
+//     D(level(x))
+//   D saves the number of polys in PS which have degree B(level(x)) in 
+//     D(level(x))
+//   E saves the minimal total degree of lcoeffs of polys wrt x in E(level(x))
+//   F saves the minimal number of terms of lcoeffs of E(level(x)) in F(~)
+//   G saves nr of poly Variable x has first deg <> 0
+
+#define __INIT_GAP__ 3
+static Varlist
+reorderb(const Varlist & difference, const CFList & PS, 
+	 const int highest_level){
+  Intarray A(1, highest_level), B(1, highest_level), C(1, highest_level), 
+    D(1, highest_level), E(1, highest_level), F(1, highest_level),
+    G(1, highest_level);
+  initArray(highest_level,A,B,C,D,E,F,G);
+  int i=0, j, n=difference.length(), gap=1+__INIT_GAP__;
+  Variable temp;
+  Array<Variable> v(0,n);
+  VarlistIterator J;
+
+  for (J= difference; J.hasItem(); J++ ){ 
+    v[i]= J.getItem(); 
+    i++ ; 
+  }
+  
+  while ( gap <= n) gap = __INIT_GAP__ * gap+1;
+  gap /= __INIT_GAP__;
+  DEBOUTLN(CERR, "gap is ", gap);
+  
+  while ( gap > 0 ){
+    for ( i=gap; i<= n-1; i++){
+      temp = v[i];
+      for ( j= i-gap; j >=0 ; j-=gap){
+	if (degord(v[j],temp, PS, A,B,C,D,E,F,G))  break;
+	v[j+gap] = v[j];
+	//CERR << "v[" << j+gap << "]= " << v[j+gap] << "\n";
+      }
+      v[j+gap] = temp;
+      //CERR << "v[" << j+gap << "]= " << v[j+gap] << "\n";
+    }
+    gap /= __INIT_GAP__;
+  }
+  Varlist output;
+  for (i=0; i<= n-1; i++)
+    output.append(v[i]);
+  DEBOUTLN(CERR, "A= ", A);  DEBOUTLN(CERR, "B= ", B);  
+  DEBOUTLN(CERR, "C= ", C);  DEBOUTLN(CERR, "D= ", D);
+  DEBOUTLN(CERR, "E= ", E);  DEBOUTLN(CERR, "F= ", F);
+  DEBOUTLN(CERR, "G= ", G);
+  return output;
+}
+
+// set up a new orderd list of Variables.
+// we try to reorder the variables heuristically optimal.
+Varlist
+neworder( const CFList & PolyList ){
+  CFList PS= PolyList, PS1=PolyList;
+  Varlist oldorder, reorder, difference;
+
+  DEBINCLEVEL(CERR, "neworder");
+  TIMING_START(neworder_time);
+
+  int highest_level= level(get_max_var(PS));
+  DEBOUTLN(CERR, "neworder: highest_level=   ", highest_level);
+  DEBOUTLN(CERR, "neworder: that is variable ", Variable(highest_level));
+
+  // set up oldorder and first criterion: only_in_one
+  for (int i=highest_level; i>=1; i--){
+    oldorder.insert(Variable(i));
+    CFList is_one= only_in_one(PS1, Variable(i)); 
+    if ( is_one.length() == 1 ){
+      DEBOUTLN(CERR, "Found a variable which is in only one Polynomial: ", 
+	       Variable(i));
+      DEBOUTLN(CERR, ".... in the Polynomial: ", is_one.getFirst());
+      reorder.insert(Variable(i));
+      PS1 = Difference(PS1, is_one);
+      DEBOUTLN(CERR, "New cutted list is: ", PS1);
+    }
+    else if ( is_one.length() == 0 ){
+      DEBOUTLN(CERR, "Found a variable which is in no polynomial: ", 
+	       Variable(i));
+      DEBOUTMSG(CERR, "... assigning it the highest level.");
+      reorder.append(Variable(i)); // assigne it the highest level
+      PS1 = Difference(PS1, is_one); 
+    }
+  }
+  difference = Difference(oldorder,reorder);
+  DEBOUTLN(CERR, "Missing variables are: ", difference);
+  // rearrange the ordering of the variables!
+  difference = reorderb(difference, PS, highest_level);
+  DEBOUTLN(CERR, "second criterion gives: ", difference);
+  reorder = Union(reorder, difference);
+  DEBOUTLN(CERR, "old order was:     ", oldorder);
+  DEBOUTLN(CERR, "New order will be: ", Union(reorder, Difference(oldorder,reorder)));
+  DEBDECLEVEL(CERR, "neworder");
+  TIMING_END(neworder_time);
+
+  TIMING_PRINT(neworder_time, "\ntime used for neworder   : ");
+
+
+  return Union(reorder,Difference(oldorder,reorder));
+}
+
+// the same as above, only returning a list of CanonicalForms
+CFList
+newordercf(const CFList & PolyList ){
+  Varlist reorder= neworder(PolyList);
+  CFList output;
+
+  for (VarlistIterator i=reorder; i.hasItem(); i++)
+    output.append(CanonicalForm(i.getItem()));
+
+  return output;
+}
+
+// the same as above, only returning a list of ints
+IntList
+neworderint(const CFList & PolyList ){
+  Varlist reorder= neworder(PolyList);
+  IntList output;
+
+  for (VarlistIterator i=reorder; i.hasItem(); i++)
+    output.append(level(i.getItem()));
+
+  return output;
+}
+
+// swapvar a whole list of CanonicalForms
+static CFList
+swapvar( const CFList & PS, const Variable & x, const Variable & y){
+  CFList ps;
+
+  for (CFListIterator i= PS; i.hasItem(); i++)
+    ps.append(swapvar(i.getItem(),x,y));
+  return ps;
+}
+
+static CFFList
+swapvar( const CFFList & PS, const Variable & x, const Variable & y){
+  CFFList ps;
+
+  for (CFFListIterator i= PS; i.hasItem(); i++)
+    ps.append(CFFactor(swapvar(i.getItem().factor(),x,y),i.getItem().exp()));
+  return ps;
+}
+
+// a library function: we reorganize the global variable ordering
+CFList
+reorder( const Varlist & betterorder, const CFList & PS){
+  int i=1, n = betterorder.length();
+  Intarray v(1,n);
+  CFList ps=PS;
+
+  //initalize:
+  for (VarlistIterator j = betterorder; j.hasItem(); j++){
+    v[i]= level(j.getItem()); i++;
+  }
+  DEBOUTLN(CERR, "reorder: Original ps=  ", ps);
+  // reorder:
+  for (i=1; i <= n; i++)
+    ps=swapvar(ps,Variable(v[i]),Variable(n+i));
+  DEBOUTLN(CERR, "reorder: Reorganized ps= ", ps); 
+  return ps;
+}
+
+CFFList
+reorder( const Varlist & betterorder, const CFFList & PS){
+  int i=1, n = betterorder.length();
+  Intarray v(1,n);
+  CFFList ps=PS;
+
+  //initalize:
+  for (VarlistIterator j = betterorder; j.hasItem(); j++){
+    v[i]= level(j.getItem()); i++;
+  }
+  DEBOUTLN(CERR, "reorder: Original ps=  ", ps);
+  // reorder:
+  for (i=1; i <= n; i++)
+    ps=swapvar(ps,Variable(v[i]),Variable(n+i));
+  DEBOUTLN(CERR, "reorder: Reorganized ps= ", ps); 
+  return ps;
+}
+
+ListCFList
+reorder(const Varlist & betterorder, const ListCFList & Q){
+  ListCFList Q1;
+
+  for ( ListCFListIterator i=Q; i.hasItem(); i++)
+    Q1.append(reorder(betterorder, i.getItem()));
+  return Q1;
+}
+
Index: factory/libfac/charset/reorder.h
===================================================================
--- factory/libfac/charset/reorder.h	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/reorder.h	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,42 @@
+////////////////////////////////////////////////////////////
+// emacs edit mode for this file is -*- C++ -*-
+////////////////////////////////////////////////////////////
+#ifndef INCL_REORDER_H
+#define INCL_REORDER_H
+
+#include <factory.h>
+#include <tmpl_inst.h>  // for typedef's
+
+/*BEGINPUBLIC*/
+// the next three give you a heuristically optimal reorderd list of the
+// variables. For internal and external (e.g. Singular/Macaulay2) library use.
+// This is really experimental!
+// See the comments in reorder.cc.
+// 
+// this gives you a heuristically optimal ordering for the ring variables
+// if you use the irreducible characteristic series.
+Varlist  neworder( const CFList & PolyList );
+// the same as neworder(...) only returning a list of CanonicalForm 's
+// (i.e. the variables as CanonicalForms)
+CFList   newordercf(const CFList & PolyList );
+// the same as neworder(...) only returning a list of int 's (i.e. the levels)
+IntList  neworderint(const CFList & PolyList );
+
+// for library internal use only:
+// next function reorders the variables in PS:
+// a code segment to use:
+// ...
+// #include <tmpl_inst.h> // for typedef's
+// CFList PS= <setup-your-list-of-CanonicalForms>;
+// Varlist betterorder= neworder(PS);
+// PS= reorder(betterorder,PS); // reorder variables in PS from oldorder 
+//                                 to betterorder
+// ListCFList Q= IrrCharSeries( PS );
+// Q= reorder(betterorder,Q);   // revert ordering to oldorder
+// 
+CFList reorder( const Varlist & betterorder, const CFList & PS);
+CFFList reorder( const Varlist & betterorder, const CFFList & PS);
+ListCFList reorder(const Varlist & betterorder, const ListCFList & Q);
+/*ENDPUBLIC*/
+
+#endif /* INCL_REORDER_H */
Index: factory/libfac/charset/tests/charsettest
===================================================================
--- factory/libfac/charset/tests/charsettest	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/charsettest	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,47 @@
+#!/bin/csh
+set OS=`uname`
+if ( $OS == "Linux" ) then
+	set ECHO="/bin/echo -e"
+	set MAKE=make
+else if ( $OS == "SunOS" ) then
+	set ECHO=/usr/5bin/echo
+	set MAKE=gmake
+else
+	set ECHO=echo
+	set MAKE=make
+endif
+
+(cd .. ; $MAKE opttestcs )
+
+$ECHO "\n\ntesting CharSet:\c"
+$ECHO "\nt1  \c" 
+../testcs < t1  
+$ECHO "\nt2  \c"
+../testcs < t2
+$ECHO "\nt3  \c"
+../testcs < t3
+$ECHO "\nt4  \c"
+../testcs < t4
+$ECHO "\nt5  \c"
+../testcs < t5
+$ECHO "\nt6  \c"
+../testcs < t6
+$ECHO "\nt7  \c"
+../testcs < t7
+$ECHO "\nt8  \c"
+../testcs < t8
+$ECHO "\nt9  \c"
+../testcs < t9
+
+$ECHO "\nw1 \c"
+../testcs < w1
+$ECHO "\nw2 \c"
+../testcs < w2
+$ECHO "\nw3 \c"
+../testcs < w3
+$ECHO "\nw4 \c"
+../testcs < w4
+
+$ECHO "\nx1 \c"
+../testcs < x1
+
Index: factory/libfac/charset/tests/hans1
===================================================================
--- factory/libfac/charset/tests/hans1	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/hans1	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,12 @@
+32003
+abcdefh
+-1*a-1*d+e+f+h;
+-2*a^2+2*a*d+2*c*e+2*b*f-1*a-1;
+3*a^3-3*a^2*d+3*c^2*e+3*b^2*f+3*a^2+4*a-1*d;
+6*a^3-6*a^2*d+6*b*e*g+6*a^2-3*a*d+4*a-1*d;
+-4*a^4+4*a^3*d+4*c^3*e+4*b^3*f-6*a^3-10*a^2+4*a*d-1*a-1;
+-8*a^4+8*a^3*d+8*b*c*e*g-12*a^3+4*a^2*d-14*a^2+4*a*d-3*a-1;
+-12*a^4+12*a^3*d+12*b^2*e*g-18*a^3+12*a^2*d-14*a^2+8*a*d-1*a-1;
+24*a^4-24*a^3*d+36*a^3-24*a^2*d+26*a^2-8*a*d+7*a+1;
+0;
+
Index: factory/libfac/charset/tests/hans2
===================================================================
--- factory/libfac/charset/tests/hans2	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/hans2	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,8 @@
+0
+abcdefghj
+d*b-e*a;
+c*e-f*b;
+g*e-h*d;
+h*f-j*e;
+0;
+
Index: factory/libfac/charset/tests/mcharsettest
===================================================================
--- factory/libfac/charset/tests/mcharsettest	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/mcharsettest	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,51 @@
+#!/bin/csh
+set OS=`uname`
+if ( $OS == "Linux" ) then
+	set ECHO="/bin/echo -e"
+	set MAKE=make
+else if ( $OS == "SunOS" ) then
+	set ECHO=/usr/5bin/echo
+	set MAKE=gmake
+else
+	set ECHO=echo
+	set MAKE=make
+endif
+
+(cd .. ; $MAKE testmcs )
+
+$ECHO "\n\ntesting CharSet:\c"
+$ECHO "\nt1  \c" 
+../testmcs < t1  
+$ECHO "\nt2  \c"
+../testmcs < t2
+$ECHO "\nt3  \c"
+../testmcs < t3
+$ECHO "\nt4  \c"
+../testmcs < t4
+$ECHO "\nt5  \c"
+../testmcs < t5
+$ECHO "\nt6  \c"
+../testmcs < t6
+$ECHO "\nt7  \c"
+../testmcs < t7
+$ECHO "\nt8  \c"
+../testmcs < t8
+$ECHO "\nt9  \c"
+../testmcs < f9
+
+$ECHO "\nw1 \c"
+../testmcs < w1
+$ECHO "\nw2 \c"
+../testmcs < w2
+$ECHO "\nw3 \c"
+../testmcs < w3
+$ECHO "\nw4 \c"
+../testmcs < w4
+
+$ECHO "\nx1 \c"
+../testmcs < x1
+
+$ECHO "\n"
+$ECHO 
+# long time section:
+$ECHO "\n Long-time section!\n"
Index: factory/libfac/charset/tests/t1
===================================================================
--- factory/libfac/charset/tests/t1	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/t1	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,8 @@
+31991
+abcd
+a*b*c-b*c*d-c*d*a-d*a*b;
+a^3-b^3-a^3*b^3+c^3+d^3;
+a*b+b*c+c*d+d*a;
+a^2;
+c*a*b+a*b+a+b;
+0;
Index: factory/libfac/charset/tests/t2
===================================================================
--- factory/libfac/charset/tests/t2	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/t2	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,6 @@
+31991
+xyz
+x^3-y^2;
+x^3-z^2;
+y^3-z^2;
+0;
Index: factory/libfac/charset/tests/t3
===================================================================
--- factory/libfac/charset/tests/t3	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/t3	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,8 @@
+31991
+abcdef
+a*b*c+a*b+a*e;
+a-b*d*e;
+d*e*f+f;
+f-e;
+a-a*b;
+0;
Index: factory/libfac/charset/tests/t4
===================================================================
--- factory/libfac/charset/tests/t4	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/t4	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,6 @@
+31991
+xyz
+x^3-x^2+y^2;
+x-y^3;
+x*y*z+z^3;
+0;
Index: factory/libfac/charset/tests/t5
===================================================================
--- factory/libfac/charset/tests/t5	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/t5	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,5 @@
+31991
+xy
+x^2+y^2-1;
+x*y;
+0;
Index: factory/libfac/charset/tests/t6
===================================================================
--- factory/libfac/charset/tests/t6	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/t6	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,10 @@
+31991
+caloqhjnbdpkmgef
+a^4*b^3*c^3-e^2*f^4+a^3*b^2*e;
+a*b*d^3-q^8;
+a*b*g^4*d^2-g^3-g^2;
+a^4-a^5-a^3+k^2*o^2*p;
+q^3*j^3*m*n;
+n^3-n^2+a*b*c;
+p^4-p^3-p^5*b^2-a*b*c*d^3-q^2*j^2;
+0;
Index: factory/libfac/charset/tests/t7
===================================================================
--- factory/libfac/charset/tests/t7	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/t7	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,6 @@
+31991
+xyzt
+t^2+x*t^2-y*t-x*y*t+x*y+3*y;
+x*t+z-x*y;
+z*t-2*y^2-x*y-1;
+0;
Index: factory/libfac/charset/tests/t8
===================================================================
--- factory/libfac/charset/tests/t8	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/t8	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,10 @@
+31991
+dcaflohjnbspkmge
+a^4*b^3*c^3-e^2*f^4+a^3*b^2*e;
+a*b*d^3-s^8;
+a*b*g^4*d^2-g^3-g^2;
+a^4-a^5-a^3+k^2*o^2*p;
+s^3*j^3*m*n;
+n^3-n^2+a*b*c;
+p^4-p^3-p^5*b^2-a*b*c*d^3-s^2*j^2;
+0;
Index: factory/libfac/charset/tests/t8a
===================================================================
--- factory/libfac/charset/tests/t8a	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/t8a	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,10 @@
+31991
+ponmlkjigfedcba
+a^4*b^3*c^3-e^2*f^4+a^3*b^2*e;
+a*b*d^3-i^8;
+a*b*g^4*d^2-g^3-g^2;
+a^4-a^5-a^3+k^2*o^2*p;
+i^3*j^3*m*n;
+n^3-n^2+a*b*c;
+p^4-p^3-p^5*b^2-a*b*c*d^3-i^2*j^2;
+0;
Index: factory/libfac/charset/tests/t9
===================================================================
--- factory/libfac/charset/tests/t9	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/t9	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,8 @@
+31991
+dcba
+a*b*c-b*c*d-c*d*a-d*a*b;
+a^3-b^3-a^3*b^3+c^3+d^3;
+a*b+b*c+c*d+d*a;
+a^2;
+c*a*b+a*b+a+b;
+0;
Index: factory/libfac/charset/tests/w1
===================================================================
--- factory/libfac/charset/tests/w1	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/w1	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,6 @@
+31991
+xyzt
+t^2+x*t^2-y*t-x*y*t+x*y+3*y;
+x*t+z-x*y;
+z*t-2*y^2-x*y-1;
+0;
Index: factory/libfac/charset/tests/w2
===================================================================
--- factory/libfac/charset/tests/w2	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/w2	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,11 @@
+31991
+acbdefgh
+h+f+e-d-a;
+2*f*b+2*e*c-1-a-2*a^2+2*a*d;
+3*f*b^2+3*e*c^2-d-3*a^2*d+4*a+3*a^2+2*a^3;
+6*e*g*b-d-3*a*d-6*a^2*d+4*a+6*a^2+6*a^3;
+4*f*b^2+4*e*c^3-1-a-10*a^2-6*a^3-4*a^4+4*d*a+4*a^3*d;
+8*e*c*g*b-1-3*a-14*a^2-12*a^3-8*a^4-4*a*d+4*a^2*d+8*a^3*d;
+12*e*c*g*b-1-a-14*a^2-18*a^3-12*a^4+8*a*d+12*a^2*d+12*a^3*d;
+1+7*a+26*a^2+36*a^3+24*a^4-8*a*d-24*a^2*d-24*a^3*d;
+0;
Index: factory/libfac/charset/tests/w3
===================================================================
--- factory/libfac/charset/tests/w3	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/w3	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,15 @@
+31991
+abcdefghkj
+-j+a^2*b^2*c^2*d^2*e+a^2*b^2*c^2*d^2*f+a^2*b^2*c^2*d*e^2+2*a^2*b^2*c^2*d*e*f+a^2*b^2*c^2*d*f^2+a^2*b^2*c^2*e^2*f+a^2*b^2*c^2*e*f^2+a^2*d*h-b^2*f*g*c^2*e*k*g-h-k;
+2*a*d*h+2*a*b^2*c^2*e*f^2+2*a*b^2*c^2*d*f^2+2*a*b^2*c^2*e^2*f+4*a*b^2*c^2*d*e*f+2*a*b^2*c^2*d^2*f+2*a*b^2*c^2*d*e^2+2*a*b^2*c^2*d^2*e;
+-2*b*c^2*e*f*g^2*k+2*a^2*b*c^2*e*f^2+2*a^2*b*c^2*d*f^2+2*a^2*b*c^2*e^2*f+4*a^2*b*c^2*d*e*f+2*a^2*b*c^2*d^2*f+2*a^2*b*c^2*d*e^2+2*a^2*b*c^2*d^2*e;
+-2*b^2*c*e*f*g^2*k+2*a^2*b^2*c*e*f^2+2*a^2*b^2*c*d*f^2+2*a^2*b^2*c*e^2*f+4*a^2*b^2*c*d*e*f+2*a^2*b^2*c*d^2*f+2*a^2*b^2*c*d*e^2+2*a^2*b^2*c*d^2*e;
+a^2*h+a^2*b^2*c^2*f^2+2*a^2*b^2*c^2*e*f+2*a^2*b^2*c^2*d*f+a^2*b^2*c^2*e^2+2*a^2*b^2*c^2*d*e;
+-b^2*c^2*f*g^2*k+a^2*b^2*c^2*f^2+2*a^2*b^2*c^2*e*f+2*a^2*b^2*c^2*d*f+2*a^2*b^2*c^2*d*e+a^2*b^2*c^2*d^2;
+-b^2*c^2*e*g^2*k+2*a^2*b^2*c^2*e*f+2*a^2*b^2*c^2*d*f+a^2*b^2*c^2*e^2+2*a^2*b^2*c^2*d*e+a^2*b^2*c^2*d^2;
+-2*b^2*c^2*e*f*g*k;
+a^2*d-1;
+-b^2*c^2*e*f*g^2-1;
+0;
+
+
Index: factory/libfac/charset/tests/w3a
===================================================================
--- factory/libfac/charset/tests/w3a	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/w3a	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,14 @@
+31991
+abcdefghij
+-6*b*c^2*e*f*g^2+2*a^2*b*c^2*e*f^2+2*a^2*b*c^2*d*f^2+2*a^2*b*c^2*e^2*f+4*a^2*b*c^2*d*e*f+2*a^2*b*c^2*d^2*f+2*a^2*b*c^2*d*e^2+2*a^2*b*c^2*d^2*e;
+-j-b^2*c^2*e*f*g^2-1+a^2*d*h-h+a^2*b^2*c^2*e*f^2+a^2*b^2*c^2*d*f^2+a^2*b^2*c^2*e^2*f+2*a^2*b^2*c^2*d*e*f+a^2*b^2*c^2*d^2*f+a^2*b^2*c^2*d*e^2+a^2*b^2*c^2*d^2*e;
+-8*b^2*c*e*f*g^2+2*a^2*b^2*c*e*f^2+2*a^2*b^2*c*d*f^2+2*a^2*b^2*c*e^2*f+4*a^2*b^2*c*d*e*f+2*a^2*b^2*c*d^2*f+2*a^2*b^2*c*d*e^2+2*a^2*b^2*c*d^2*e;
+-16*b^2*c^2*e*f*g;
+-7*b^2*c^2*e*g^2+2*a^2*b^2*c^2*e*f+2*a^2*b^2*c^2*d*f+a^2*b^2*c^2*e^2+2*a^2*b^2*c^2*d*e+a^2*b^2*c^2*d^2;
+-6*b^2*c^2*f*g^2+a^2*b^2*c^2*f^2+2*a^2*b^2*c^2*e*f+2*a^2*b^2*c^2*d*f+2*a^2*b^2*c^2*d*e+a^2*b^2*c^2*d^2;
+a^2*d-1;
+-b^2*c^2*e*f*g^2-1;
+2*a*d*h+2*a*b^2*c^2*e*f^2+2*a*b^2*c^2*d*f^2+2*a*b^2*c^2*e^2*f+4*a*b^2*c^2*d*e*f+2*a*b^2*c^2*d^2*f+2*a*b^2*c^2*d*e^2+2*a*b^2*c^2*d^2*e;
+a^2*h+a^2*b^2*c^2*f^2+2*a^2*b^2*c^2*e*f+2*a^2*b^2*c^2*d*f+a^2*b^2*c^2*e^2+2*a^2*b^2*c^2*d*e;
+0;
+
Index: factory/libfac/charset/tests/w4
===================================================================
--- factory/libfac/charset/tests/w4	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/w4	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,6 @@
+31991
+zyxr
+x^2+y^2+z^2-r^2;
+x*y+z^2-1;
+x*y*z-x^2-y^2-z+1;
+0;
Index: factory/libfac/charset/tests/x1
===================================================================
--- factory/libfac/charset/tests/x1	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/x1	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,6 @@
+31991
+abcxyz
+a*(x+y+z)-(x-y-z);
+b*(x+z)-(x+y-z);
+c*(x-y+z)-(x+y);
+0;
Index: factory/libfac/charset/tests/x2
===================================================================
--- factory/libfac/charset/tests/x2	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/x2	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,23 @@
+31991
+abcdefghjklmnopqr
+d*k;
+e*g+d*j+c*k+l;
+b*h;
+e*l;
+b*f+b*g+a*h+d*h+b*k+h+o;
+a*f+d*f+a*g+d*g+c*h+e*h+b*j+a*k+d*k+b*l+f+g+k+m+n+q;
+c*f+e*f+c*g+e*g+a*j+d*j+c*k+e*k+a*l+d*l+j+l+p+r-1;
+c*j+e*j+c*l+e*l;
+e*j+c*l+2*e*l;
+e*f+e*g+d*j+c*k+2*e*k+a*l+2*d*l+l+r;
+d*f+d*g+e*h+a*k+2*d*k+b*l+k+q;
+d*h+b*k;
+e*k+d*l;
+2*c*j+e*j+c*l;
+c*f+2*c*g+e*g+a*j+d*j+c*k+2*j+l+p;
+a*g+d*g+c*h+b*j+f+2*g+k+n;
+b*g+h;
+e*j+c*l;
+d*g+k;
+0;
+
Index: factory/libfac/charset/tests/x3
===================================================================
--- factory/libfac/charset/tests/x3	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/x3	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,13 @@
+31991
+rstuvwxyz
+s*u+r*v;
+t*u+r*w;
+t*v+s*w;
+s*x+r*y;
+t*x+r*z;
+t*y+s*z;
+v*x+u*y;
+w*x+u*z;
+w*y+v*z;
+0;
+
Index: factory/libfac/charset/tests/x4
===================================================================
--- factory/libfac/charset/tests/x4	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/x4	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,36 @@
+31991
+abcdef
+a*d*e*f+3/2*b*e^2*f-1/2*c*e*f^2;
+    a*d^2*f+5/4*b*d*e*f+1/4*c*d*f^2;
+    a*d^2*e+3/4*b*d*e^2+7/4*c*d*e*f;
+    a*c*d*e+3/2*b*c*e^2-1/2*c^2*e*f;
+    a*c*d^2+5/4*b*c*d*e+1/4*c^2*d*f;
+    a*b*d*f+3/2*b^2*e*f-1/2*b*c*f^2;
+    a*b*d*e+3/4*b^2*e^2-a*c*d*f+1/2*b*c*e*f-1/4*c^2*f^2;
+    a*b*d^2+3/4*b^2*d*e+7/4*b*c*d*f;
+    a*b*c*d+3/2*b^2*c*e-1/2*b*c^2*f;
+    a^2*d*f+5/4*a*b*e*f+1/4*a*c*f^2;
+    a^2*d*e+3/4*a*b*e^2+7/4*a*c*e*f;
+    a^2*d^2-9/16*b^2*e^2+2*a*c*d*f-9/8*b*c*e*f+7/16*c^2*f^2;
+    a^2*c*d+5/4*a*b*c*e+1/4*a*c^2*f;
+    a^2*b*d+3/4*a*b^2*e+7/4*a*b*c*f;
+    b*c^3*d+1/4*c*d^3*e;
+    b^2*c^2*e-1/3*b*c^3*f+2/3*c*d^2*e*f;
+    b^2*c^2*d-1/2*a*d^4-3/8*b*d^3*e-1/8*c*d^3*f;
+    b^3*c*e-1/3*b^2*c^2*f+2/3*b*d^2*e*f;
+    b^3*c*d+3/4*b*d^3*f;
+    a*c^3*e-1/3*c*d*e^3;
+    a*c^3*d-3/4*b*c^3*e-1/2*c*d^2*e^2+1/4*c^4*f;
+    a*b*c^2*f-c*d*e*f^2;
+    a*b*c^2*e-c*d*e^2*f;
+    a*b^2*c*f-b*d*e*f^2;
+    a*b^3*f-3*b*d*f^3;
+    a*b^3*d+3/4*b^4*e-1/4*b^3*c*f-3/2*b*d^2*f^2;
+    a^2*c^2*e-2/3*a*d*e^3-1/2*b*e^4+1/6*c*e^3*f;
+    a^2*b*c*f+3/2*b*e^2*f^2-1/2*c*e*f^3;
+    a^2*b*c*e+3/2*b*e^3*f-1/2*c*e^2*f^2;
+    a^2*b^2*f-2*a*d*f^3+3/2*b*e*f^3-1/2*c*f^4;
+    a^3*c*e+4/3*a*e^3*f;
+    a^3*b*f+4*a*e*f^3;
+    a^4*d+3/4*a^3*b*e+1/4*a^3*c*f-2*a*e^2*f^2;
+0;
Index: factory/libfac/charset/tests/x5
===================================================================
--- factory/libfac/charset/tests/x5	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/x5	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,11 @@
+32003
+abcdefgh
+h+f+e-(d+a);
+2*f*b+2*e*c-(1+a+2*a^2-2*d*a);
+3*f*b^2+3*e*c^2-(d*(1+3*a^2)-4*a-3*a^2-3*a^3);
+6*e*g*b-(d*(1+3*a+6*a^2)-4*a-6*a^2-6*a^3);
+4*f*b^3+4*e*c^3-(1+a+10*a^2+6*a^3+4*a^4-4*d*(a+a^3));
+8*e*c*g*b-(1+3*a+14*a^2+12*a^3+8*a^4-d*(4*a+4*a^2+8*a^3));
+12*e*g*b^2-(1+a+14*a^2+18*a^3+12*a^4-d*(8*a+12*a^2+12*a^3));
+1+7*a+26*a^2+36*a^3+24*a^4-d*(8*a+24*a^2+24*a^3);
+0;
Index: factory/libfac/charset/tests/x9
===================================================================
--- factory/libfac/charset/tests/x9	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
+++ factory/libfac/charset/tests/x9	(revision 16055bd66cbe17a643cdc139b5ad13307053b068)
@@ -0,0 +1,6 @@
+31991
+xy
+x^2;
+x*y;
+0;
+
