Index: factory/facAlgFunc.cc
===================================================================
--- factory/facAlgFunc.cc	(revision 24819474a270d2480ed2d9f37a85ab3f5eafee1c)
+++ factory/facAlgFunc.cc	(revision 5355082c234c56e9e117c19600ae46cf6324f5ca)
@@ -3,26 +3,142 @@
 ////////////////////////////////////////////////////////////
 
-// FACTORY - Includes
-#include <factory.h>
-// Factor - Includes
-#include <tmpl_inst.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"
+#include "config.h"
+
+#include "cf_assert.h"
+#include "debug.h"
+
+#include "algext.h"
+#include "cf_random.h"
+#include "cf_generator.h"
+#include "cf_irred.h"
+#include "cf_iter.h"
+#include "cf_util.h"
+#include "cf_algorithm.h"
+#include "cf_map.h"
+#include "cfModResultant.h"
+#include "cfCharSets.h"
+#include "facAlgFunc.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
+CanonicalForm
+Prem (const CanonicalForm& F, const CanonicalForm& G)
+{
+  CanonicalForm f, g, l, test, lu, lv, t, retvalue;
+  int degF, degG, levelF, levelG;
+  bool reord;
+  Variable v, vg= G.mvar();
+
+  if ( (levelF= F.level()) < (levelG= G.level()))
+    return F;
+  else
+  {
+    if ( levelF == levelG )
+    {
+      f= F;
+      g= G;
+      reord= false;
+      v= F.mvar();
+    }
+    else
+    {
+      v= Variable (levelF + 1);
+      f= swapvar (F, vg, v);
+      g= swapvar (G, vg, v);
+      reord= true;
+    }
+    degG= degree (g, v );
+    degF= degree (f, v );
+    if (degG <= degF)
+    {
+      l= LC (g);
+      g= g - l*power (v, degG);
+    }
+    else
+      l= 1;
+    while ((degG <= degF) && (!f.isZero()))
+    {
+      test= gcd (l, LC(f));
+      lu= l / test;
+      lv= LC(f) / test;
+
+      t= power (v, degF - degG)*g*lv;
+
+      if (degF == 0)
+        f= 0;
+      else
+        f= f - LC (f)*power (v, degF);
+
+      f= lu*f - t;
+      degF= degree (f, v);
+    }
+
+    if (reord)
+      retvalue= swapvar (f, vg, v);
+    else
+      retvalue= f;
+
+    return retvalue;
+  }
+}
+
+static
+CanonicalForm
+Prem( const CanonicalForm &f, const CFList &L)
+{
+  CanonicalForm rem= f;
+  CFListIterator i= L;
+
+  for (i.lastItem(); i.hasItem(); i--)
+    rem= normalize (Prem (rem, i.getItem()));
+
+  return rem;
+}
+
+static
+CanonicalForm
+Prem (const CanonicalForm &f, const CFList &L, const CFList &as)
+{
+  CanonicalForm rem = f;
+  CFListIterator i = L;
+  for ( i.lastItem(); i.hasItem(); i-- )
+    rem = Prem( rem, i.getItem() );
+  return normalize (rem); //TODO better normalize in case as.length() == 1 && as.getFirst().level() == 1 ???
+}
+
+CFFList
+myappend( const CFFList & Inputlist, const CFFactor & TheFactor)
+{
+  CFFList Outputlist ;
+  CFFactor copy;
+  CFFListIterator i;
+  int exp=0;
+
+  for ( i=Inputlist ; i.hasItem() ; i++ )
+  {
+    copy = i.getItem();
+    if ( copy.factor() == TheFactor.factor() )
+      exp += copy.exp();
+    else
+      Outputlist.append(copy);
+  }
+  Outputlist.append( CFFactor(TheFactor.factor(), exp + TheFactor.exp()));
+  return Outputlist;
+}
+
+CFFList
+myUnion(const CFFList & Inputlist1,const CFFList & Inputlist2)
+{
+  CFFList Outputlist;
+  CFFListIterator i;
+
+  for ( i=Inputlist1 ; i.hasItem() ; i++ )
+    Outputlist = myappend(Outputlist, i.getItem() );
+  for ( i=Inputlist2 ; i.hasItem() ; i++ )
+    Outputlist = myappend(Outputlist, i.getItem() );
+
+  return Outputlist;
+}
 
 ///////////////////////////////////////////////////////////////
@@ -34,6 +150,6 @@
 {
   FFRandom gen;
-  if (degree (Extension) < 0)
-    factoryError("libfac: evaluate: Extension not inFF() or inGF() !");
+  /*if (degree (Extension) < 0)
+    factoryError("libfac: evaluate: Extension not inFF() or inGF() !");*/
   return find_irreducible( degree_of_Extension, gen, Variable(1) );
 }
@@ -64,4 +180,321 @@
 {
     current++;
+}
+
+int hasAlgVar(const CanonicalForm &f, const Variable &v)
+{
+  if (f.inBaseDomain()) return 0;
+  if (f.inCoeffDomain())
+  {
+    if (f.mvar()==v) return 1;
+    return hasAlgVar(f.LC(),v);
+  }
+  if (f.inPolyDomain())
+  {
+    if (hasAlgVar(f.LC(),v)) return 1;
+    for( CFIterator i=f; i.hasTerms(); i++)
+    {
+      if (hasAlgVar(i.coeff(),v)) return 1;
+    }
+  }
+  return 0;
+}
+
+int hasVar(const CanonicalForm &f, const Variable &v)
+{
+  if (f.inBaseDomain()) return 0;
+  if (f.inCoeffDomain())
+  {
+    if (f.mvar()==v) return 1;
+    return hasAlgVar(f.LC(),v);
+  }
+  if (f.inPolyDomain())
+  {
+    if (f.mvar()==v) return 1;
+    if (hasVar(f.LC(),v)) return 1;
+    for( CFIterator i=f; i.hasTerms(); i++)
+    {
+      if (hasVar(i.coeff(),v)) return 1;
+    }
+  }
+  return 0;
+}
+
+int hasAlgVar(const CanonicalForm &f)
+{
+  if (f.inBaseDomain()) return 0;
+  if (f.inCoeffDomain())
+  {
+    if (f.level()!=0)
+      return 1;
+    return hasAlgVar(f.LC());
+  }
+  if (f.inPolyDomain())
+  {
+    if (hasAlgVar(f.LC())) return 1;
+    for( CFIterator i=f; i.hasTerms(); i++)
+    {
+      if (hasAlgVar(i.coeff())) return 1;
+    }
+  }
+  return 0;
+}
+
+/// pseudo division of f and g wrt. x s.t. multiplier*f=q*g+r
+void
+psqr (const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, 
+      CanonicalForm & r, CanonicalForm& multiplier, const Variable& x)
+{
+    ASSERT( x.level() > 0, "type error: polynomial variable expected" );
+    ASSERT( ! g.isZero(), "math error: division by zero" );
+
+    // swap variables such that x's level is larger or equal
+    // than both f's and g's levels.
+    Variable X;
+    if (f.level() > g.level())
+      X= f.mvar();
+    else
+      X= g.mvar();
+    if (X.level() < x.level())
+      X= x;
+    CanonicalForm F= swapvar (f, x, X);
+    CanonicalForm G= swapvar (g, x, X);
+
+    // now, we have to calculate the pseudo remainder of F and G
+    // w.r.t. X
+    int fDegree= degree (F, X);
+    int gDegree= degree (G, X);
+    if (fDegree < 0 || fDegree < gDegree)
+    {
+      q= 0;
+      r= f;
+    }
+    else
+    {
+      CanonicalForm LCG= LC (G, X);
+      multiplier= power (LCG, fDegree - gDegree + 1);
+      divrem (multiplier*F, G, q, r);
+      q= swapvar (q, x, X);
+      r= swapvar (r, x, X);
+    }
+}
+
+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;
+
+  if (f.inCoeffDomain())
+  {
+    bool isRat=isOn(SW_RATIONAL);
+    if (getCharacteristic() == 0)
+      On(SW_RATIONAL);
+    q=ff/f;
+    if (!isRat && getCharacteristic() == 0)
+      Off(SW_RATIONAL);
+  }
+  else
+    r= Sprem(ff,f,m,q);
+
+  r= Prem(q,as);
+  return r;
+}
+
+CanonicalForm
+alg_content (const CanonicalForm& f, const CFList& as)
+{
+  if (!f.inCoeffDomain())
+  {
+    CFIterator i= f;
+    CanonicalForm result= abs (i.coeff());
+    i++;
+    while (i.hasTerms() && !result.isOne())
+    {
+      result= alg_gcd (i.coeff(), result, as);
+      i++;
+    }
+    return result;
+  }
+
+  return abs (f);
+}
+
+CanonicalForm alg_gcd(const CanonicalForm & fff, const CanonicalForm &ggg,
+                      const CFList &as)
+{
+  if (fff.inCoeffDomain() || ggg.inCoeffDomain())
+    return 1;
+  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;
+  }
+
+  int v= as.getLast().level();
+  if (f.level() <= v || g.level() <= v)
+    return 1;
+
+  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;
+    if (hasVar (g, v))
+      has_alg_var=true;
+  }
+  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*/
+    if (!hasAlgVar (f) && !hasAlgVar (g))
+      res=gcd(f,g);
+
+    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())
+    return CanonicalForm(1);
+
+  CanonicalForm c_f= alg_content (f, as);
+
+  if (mvf != mvg)
+  {
+    res= alg_gcd (g, c_f, as);
+    return res;
+  }
+  Variable x= f.mvar();
+
+  // now: mvf==mvg, f.level()==g.level()
+  CanonicalForm c_g= alg_content (g, as);
+
+  int delta= degree (f) - degree (g);
+
+  f= divide (f, c_f, as);
+  g= divide (g, c_g, as);
+
+  // gcd of contents
+  CanonicalForm c_gcd= alg_gcd (c_f, c_g, as);
+  CanonicalForm tmp;
+
+  if (delta < 0)
+  {
+    tmp= f;
+    f= g;
+    g= tmp;
+    delta= -delta;
+  }
+
+  CanonicalForm r= 1;
+
+  while (degree (g, x) > 0)
+  {
+    r= Prem (f, g, as);
+    r= Prem (r, as);
+    if (!r.isZero())
+    {
+      r= divide (r, alg_content (r,as), as);
+      r /= vcontent (r,Variable (v+1));
+    }
+    f= g;
+    g= r;
+  }
+
+  if (degree (g, x) == 0)
+    return c_gcd;
+
+  c_f= alg_content (f, as);
+
+  f= divide (f, c_f, as);
+
+  f *= c_gcd;
+  f /= vcontent (f, Variable (v+1));
+
+  return f;
 }
 
@@ -112,16 +545,11 @@
   CFFListIterator i;
 
-  DEBOUTLN(CERR, "sqrf_norm_sub:      f= ", f);
-  DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);
   myrandom.reset();   s=myrandom.item();   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);
     R /= content (R);
     // sqfree check ; R is a polynomial in K[x]
@@ -129,18 +557,11 @@
     {
       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 (hasFirstAlgVar(R,X))
@@ -148,5 +569,5 @@
       else
         testlist= factorize(R);
-      DEBOUTLN(CERR, "testlist= ", testlist);
+
       if (testlist.getFirst().factor().inCoeffDomain())
         testlist.removeFirst();
@@ -160,10 +581,8 @@
         }
       }
-      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()));
@@ -171,7 +590,5 @@
         t= CanonicalForm(myrandom.item());
       s= t;
-      DEBOUTLN(CERR, "sqrf_norm_sub: testing s= ", s);
       g= f(f.mvar()-t*Palpha.mvar(), f.mvar());
-      DEBOUTLN(CERR, "             gives g= ", g);
     }
   }
@@ -188,16 +605,11 @@
   CFFListIterator i;
 
-  DEBOUTLN(CERR, "sqrf_norm_sub:      f= ", f);
-  DEBOUTLN(CERR, "sqrf_norm_sub: Palpha= ", Palpha);
   myrandom.reset();   s=myrandom.item();   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);
     R /= content (R);
     // sqfree check ; R is a polynomial in K[x]
@@ -205,18 +617,11 @@
     {
       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 (hasFirstAlgVar(R,X))
@@ -224,5 +629,4 @@
       else
         testlist= factorize(R);
-      DEBOUTLN(CERR, "testlist= ", testlist);
       if (testlist.getFirst().factor().inCoeffDomain())
         testlist.removeFirst();
@@ -236,10 +640,8 @@
         }
       }
-      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()));
@@ -247,7 +649,5 @@
         t= CanonicalForm(myrandom.item());
       s= t;
-      DEBOUTLN(CERR, "sqrf_norm_sub: testing s= ", s);
       g= f(f.mvar()-t*Palpha.mvar(), f.mvar());
-      DEBOUTLN(CERR, "             gives g= ", g);
     }
   }
@@ -259,20 +659,11 @@
            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);
@@ -281,5 +672,4 @@
   {
     FFGenerator myrandom;
-    DEBOUTMSG(CERR, "sqrf_norm: degree of extension is 0");
     sqrf_norm_sub(f,PPalpha, myrandom, s,g,R);
   }
@@ -347,101 +737,11 @@
       l= l+1;
       if ( igcd(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; }
+      else { break; }
     }
     k= k+1; l=0;
   }
   while ( 1 );
-}
-
-
-/// pseudo division of f and g wrt. x s.t. multiplier*f=q*g+r
-/// but only if the leading coefficient of g is of level lower than coeffLevel
-void
-psqr (const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q,
-      CanonicalForm & r, CanonicalForm& multiplier, const Variable& x,
-      int coeffLevel)
-{
-  ASSERT( x.level() > 0, "type error: polynomial variable expected" );
-  ASSERT( ! g.isZero(), "math error: division by zero" );
-
-  // swap variables such that x's level is larger or equal
-  // than both f's and g's levels.
-  Variable X;
-  if (f.level() > g.level())
-    X= f.mvar();
-  else
-    X= g.mvar();
-  if (X.level() < x.level())
-    X= x;
-  CanonicalForm F= swapvar (f, x, X);
-  CanonicalForm G= swapvar (g, x, X);
-
-  // now, we have to calculate the pseudo remainder of F and G
-  // w.r.t. X
-  int fDegree= degree (F, X);
-  int gDegree= degree (G, X);
-  if (fDegree < 0 || fDegree < gDegree)
-  {
-    q= 0;
-    r= f;
-  }
-  else
-  {
-    CanonicalForm LCG= LC (G, X);
-    if (LCG.level() < coeffLevel)
-    {
-      multiplier= power (LCG, fDegree - gDegree + 1);
-      divrem (multiplier*F, G, q, r);
-      q= swapvar (q, x, X);
-      r= swapvar (r, x, X);
-    }
-    else
-    {
-      q= 0;
-      r= f;
-    }
-  }
-}
-
-/// pseudo division of f and g wrt. x s.t. multiplier*f=q*g+r
-void
-psqr (const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, 
-      CanonicalForm & r, CanonicalForm& multiplier, const Variable& x)
-{
-    ASSERT( x.level() > 0, "type error: polynomial variable expected" );
-    ASSERT( ! g.isZero(), "math error: division by zero" );
-
-    // swap variables such that x's level is larger or equal
-    // than both f's and g's levels.
-    Variable X;
-    if (f.level() > g.level())
-      X= f.mvar();
-    else
-      X= g.mvar();
-    if (X.level() < x.level())
-      X= x;
-    CanonicalForm F= swapvar (f, x, X);
-    CanonicalForm G= swapvar (g, x, X);
-
-    // now, we have to calculate the pseudo remainder of F and G
-    // w.r.t. X
-    int fDegree= degree (F, X);
-    int gDegree= degree (G, X);
-    if (fDegree < 0 || fDegree < gDegree)
-    {
-      q= 0;
-      r= f;
-    }
-    else
-    {
-      CanonicalForm LCG= LC (G, X);
-      multiplier= power (LCG, fDegree - gDegree + 1);
-      divrem (multiplier*F, G, q, r);
-      q= swapvar (q, x, X);
-      r= swapvar (r, x, X);
-    }
 }
 
@@ -568,7 +868,4 @@
   bool isRat= isOn (SW_RATIONAL);
 
-  DEBOUTLN(CERR, "simpleextension: Astar= ", Astar);
-  DEBOUTLN(CERR, "simpleextension:     R= ", R);
-  DEBOUTLN(CERR, "simpleextension: Extension= ", Extension);
   CFListIterator j;
   if (Astar.length() == 1)
@@ -654,7 +951,4 @@
       }
 
-      DEBOUTLN(CERR, "simpleextension: g= ", g);
-      DEBOUTLN(CERR, "simpleextension: s= ", s);
-      DEBOUTLN(CERR, "simpleextension: R= ", R);
       Returnlist.append (ra);
       if (isFunctionField)
@@ -680,5 +974,5 @@
     return alg_lc(f.LC());
   }
-  //assert(f.inCoeffDomain());
+
   return f;
 }
@@ -851,11 +1145,6 @@
 
 // the heart of the algorithm: the one from Trager
-#ifndef DEBUGOUTPUT
-static CFFList
-alg_factor( const CanonicalForm & F, const CFList & Astar, const Variable & vminpoly, const Varlist /*& oldord*/, const CFList & as, bool isFunctionField)
-#else
 static CFFList
 alg_factor( const CanonicalForm & F, const CFList & Astar, const Variable & vminpoly, const Varlist & oldord, const CFList & as, bool isFunctionField)
-#endif
 {
   bool isRat= isOn (SW_RATIONAL);
@@ -864,11 +1153,5 @@
   CFList substlist, backSubsts;
 
-  DEBINCLEVEL(CERR,"alg_factor");
-  DEBOUTLN(CERR, "alg_factor: f= ", f);
-
   substlist= simpleextension (backSubsts, Astar, vminpoly, isFunctionField, Rstar);
-  DEBOUTLN(CERR, "alg_factor: substlist= ", substlist);
-  DEBOUTLN(CERR, "alg_factor: minpoly Rstar= ", Rstar);
-  DEBOUTLN(CERR, "alg_factor: vminpoly= ", vminpoly);
 
   f= subst (f, Astar, substlist, Rstar, isFunctionField);
@@ -902,5 +1185,4 @@
   // after here we are over an extension of a function field
 
-
   // make quasi monic
   CFList Rstarlist= CFList (Rstar);
@@ -927,16 +1209,9 @@
 
   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);
+
   Variable X;
   if (hasFirstAlgVar(R,X))
   {
     // factorize R over alg.extension with X
-    //CERR << "alg: "<< X << " mipo=" << getMipo(X,Variable('X')) <<"\n";
-    DEBOUTLN(CERR, "alg_factor: factorize( ", R);
     Factorlist =  factorize( R, X );
   }
@@ -944,9 +1219,7 @@
   {
     // factor R over k
-    DEBOUTLN(CERR, "alg_factor: Factorize( ", R);
     Factorlist = factorize(R);
   }
 
-  DEBOUTLN(CERR, "alg_factor: Factorize(R)= ", Factorlist);
   if ( !Factorlist.getFirst().factor().inCoeffDomain() )
     Factorlist.insert(CFFactor(1,1));
@@ -963,5 +1236,4 @@
   {
     g= f;
-    DEBOUTLN(CERR, "alg_factor: g= ", g);
     for ( CFFListIterator i=Factorlist; i.hasItem(); i++)
     {
@@ -974,6 +1246,4 @@
         fnew= reduce (fnew, Rstar);
       }
-
-      DEBOUTLN(CERR, "alg_factor: fnew= ", fnew);
 
       h= alg_gcd (g, fnew, Rstarlist);
@@ -1018,10 +1288,5 @@
     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");
+
   if (!isRat && getCharacteristic() == 0)
     Off (SW_RATIONAL);
@@ -1225,7 +1490,6 @@
 
   // factor F with minimal polys given in asnew
-  PremForm rem;
   asnew.append (F);
-  asnew= mcharset (asnew, rem); // TODO use modCharSet
+  asnew= charSetViaModCharSet (asnew); // TODO use modCharSet
 
   F= asnew.getLast();
@@ -1269,10 +1533,9 @@
   for (CFFListIterator k= tmp; k.hasItem(); k++)
   {
-    rem= PremForm();
     transform= transBack;
     CanonicalForm factor= k.getItem().factor();
     factor= M (factor);
     transform.append (factor);
-    transform= mcharset (transform, rem); //TODO use modCharSet
+    transform= charSetViaModCharSet (transform); //TODO use modCharSet
     for (i= transform; i.hasItem(); i++)
     {
@@ -1341,18 +1604,10 @@
 
   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()))
+  if (vf.level() <= as.getLast().level())
   {
 // ||( (as.length()==1) && (degree(f,vf)==3) && (degree(as.getFirst()==2)) )
-    DEBDECLEVEL(CERR,"newfactoras");
     if (!isRat && getCharacteristic() == 0)
       Off (SW_RATIONAL);
@@ -1379,11 +1634,7 @@
   }
   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 proved irr. of f over no extensions
   if ( Astar.length() == 0 ){
-    DEBDECLEVEL(CERR,"newfactoras");
     if (!isRat && getCharacteristic() == 0)
       Off (SW_RATIONAL);
@@ -1399,5 +1650,4 @@
 //   polynomials, we don't have transzendentals.
   Varlist newuord=Var_is_in_AS(uord,Astar);
-  DEBOUTLN(CERR, "newuord is: ", newuord);
 
   CFFList Factorlist;
@@ -1410,10 +1660,7 @@
     Fgcd= alg_gcd(f,f.deriv(),Astar);
 
-  if ( Fgcd == 0 ) {DEBOUTMSG(CERR, "WARNING: p'th root ?");}
-
   bool derivZero= f.deriv().isZero();
   if (isFunctionField && (degree (Fgcd, f.mvar()) > 0) && !derivZero)
   {
-    DEBOUTLN(CERR, "Nontrivial GCD found of ", f);
     CanonicalForm Ggcd= divide(f, Fgcd,Astar);
     if (getCharacteristic() == 0)
@@ -1425,8 +1672,6 @@
       return result;
     }
-    DEBOUTLN(CERR, "  split into ", Fgcd);
-    DEBOUTLN(CERR, "         and ", Ggcd);
+
     Fgcd= pp(Fgcd); Ggcd= pp(Ggcd);
-    DEBDECLEVEL(CERR,"newfactoras");
     if (!isRat && getCharacteristic() == 0)
       Off (SW_RATIONAL);
@@ -1441,45 +1686,32 @@
     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, isFunctionField);
-      DEBDECLEVEL(CERR,"newfactoras");
       return Factorlist;
     }
     else if ( inseperable(Astar) > 0 || derivZero) // Look if extensions are seperable
     {
-      DEBOUTMSG(CERR, "Inseperable extensions! Using Endler!");
       Factorlist= SteelTrager(f, Astar, newuord);
-      DEBOUTLN(CERR, "Endler gives: ", Factorlist);
       return Factorlist;
     }
     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, isFunctionField);
-      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, isFunctionField);
-    DEBDECLEVEL(CERR,"newfactoras");
     if (!isRat && getCharacteristic() == 0)
       Off (SW_RATIONAL);
@@ -1487,5 +1719,4 @@
   }
 
-  DEBDECLEVEL(CERR,"newfactoras");
   return CFFList(CFFactor(f,1));
 }
@@ -1508,5 +1739,5 @@
     return Factors;
   }
-  if ( cls(f) <= cls(as.getLast()) )
+  if (f.level() <= as.getLast().level())
   {
     if (!isRat && getCharacteristic() == 0)
Index: factory/facAlgFunc.h
===================================================================
--- factory/facAlgFunc.h	(revision 24819474a270d2480ed2d9f37a85ab3f5eafee1c)
+++ factory/facAlgFunc.h	(revision 5355082c234c56e9e117c19600ae46cf6324f5ca)
@@ -4,9 +4,8 @@
 ////////////////////////////////////////////////////////////
 
-#ifndef INCL_NEW_ALGFACTOR_H
-#define INCL_NEW_ALGFACTOR_H
+#ifndef FAC_ALG_FUNC_H
+#define FAC_ALG_FUNC_H
 
-#include <factory.h>
-#include <tmpl_inst.h>  // for typedef's
+#include "canonicalform.h"
 
 // missing class: IntGenerator:
@@ -26,6 +25,7 @@
 };
 
+CanonicalForm alg_gcd(const CanonicalForm &, const CanonicalForm &, const CFList &);
+/*BEGINPUBLIC*/
 CFFList newfactoras( const CanonicalForm & f, const CFList & as, int &success);
-/*BEGINPUBLIC*/
 CFFList newcfactor(const CanonicalForm & f, const CFList & as, int & success );
 /*ENDPUBLIC*/
Index: factory/libfac/charset/csutil.cc
===================================================================
--- factory/libfac/charset/csutil.cc	(revision 24819474a270d2480ed2d9f37a85ab3f5eafee1c)
+++ factory/libfac/charset/csutil.cc	(revision 5355082c234c56e9e117c19600ae46cf6324f5ca)
@@ -11,67 +11,4 @@
 extern void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
 extern CanonicalForm alg_lc(const CanonicalForm &f);
-
-int hasAlgVar(const CanonicalForm &f, const Variable &v)
-{
-  if (f.inBaseDomain()) return 0;
-  if (f.inCoeffDomain())
-  {
-    if (f.mvar()==v) return 1;
-    return hasAlgVar(f.LC(),v);
-  }
-  if (f.inPolyDomain())
-  {
-    if (hasAlgVar(f.LC(),v)) return 1;
-    for( CFIterator i=f; i.hasTerms(); i++)
-    {
-      if (hasAlgVar(i.coeff(),v)) return 1;
-    }
-  }
-  return 0;
-}
-
-int hasVar(const CanonicalForm &f, const Variable &v)
-{
-  if (f.inBaseDomain()) return 0;
-  if (f.inCoeffDomain())
-  {
-    if (f.mvar()==v) return 1;
-    return hasAlgVar(f.LC(),v);
-  }
-  if (f.inPolyDomain())
-  {
-    if (f.mvar()==v) return 1;
-    if (hasVar(f.LC(),v)) return 1;
-    for( CFIterator i=f; i.hasTerms(); i++)
-    {
-      if (hasVar(i.coeff(),v)) return 1;
-    }
-  }
-  return 0;
-}
-
-int hasAlgVar(const CanonicalForm &f)
-{
-  if (f.inBaseDomain()) return 0;
-  if (f.inCoeffDomain())
-  {
-    if (f.level()!=0)
-    {
-      //CERR << "hasAlgVar:" << f.mvar() <<"\n";
-      return 1;
-    }
-    return hasAlgVar(f.LC());
-  }
-  if (f.inPolyDomain())
-  {
-    if (hasAlgVar(f.LC())) return 1;
-    for( CFIterator i=f; i.hasTerms(); i++)
-    {
-      if (hasAlgVar(i.coeff())) return 1;
-    }
-  }
-  return 0;
-}
-
 
 static bool
@@ -278,29 +215,4 @@
 }
 
-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 isRat=isOn(SW_RATIONAL);
-    if (getCharacteristic() == 0)
-      On(SW_RATIONAL);
-    q=ff/f;
-    if (!isRat && getCharacteristic() == 0)
-      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).
@@ -870,151 +782,4 @@
 #endif
 
-CanonicalForm
-alg_content (const CanonicalForm& f, const CFList& as)
-{
-  if (!f.inCoeffDomain())
-  {
-    CFIterator i= f;
-    CanonicalForm result= abs (i.coeff());
-    i++;
-    while (i.hasTerms() && !result.isOne())
-    {
-      result= alg_gcd (i.coeff(), result, as);
-      i++;
-    }
-    return result;
-  }
-
-  return abs (f);
-}
-
-CanonicalForm alg_gcd(const CanonicalForm & fff, const CanonicalForm &ggg,
-                      const CFList &as)
-{
-  if (fff.inCoeffDomain() || ggg.inCoeffDomain())
-    return 1;
-  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(",fff," \n ");
-  //out_cf("ggg",ggg,")\n");
-  int v= as.getLast().level();
-  if (f.level() <= v || g.level() <= v)
-    return 1;
-
-  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:
-  CanonicalForm c_f= alg_content (f, as);
-
-  //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_f, as);
-    return res;
-  }
-  Variable x= f.mvar();
-  // now: mvf==mvg, f.level()==g.level()
-  // content of g
-  CanonicalForm c_g= alg_content (g, as);
-
-  int delta= degree (f) - degree (g);
-
-  //f/=c_gcd;
-  //g/=c_gcd;
-  f= divide (f, c_f, as);
-  g= divide (g, c_g, as);
-
-  // gcd of contents
-  CanonicalForm c_gcd= alg_gcd (c_f, c_g, as);
-  CanonicalForm tmp;
-
-  if (delta < 0)
-  {
-    tmp= f;
-    f= g;
-    g= tmp;
-    delta= -delta;
-  }
-
-  CanonicalForm r=1;
-
-  while (degree (g, x) > 0)
-  {
-    r= Prem (f, g, as);
-    r= Prem (r, as);
-    if (!r.isZero())
-    {
-      r= divide (r, alg_content (r,as), as);
-      r /= vcontent (r,Variable (v+1));
-    }
-    f= g;
-    g= r;
-  }
-
-  if (degree (g, x) == 0)
-    return c_gcd;
-
-  c_f= alg_content (f, as);
-
-  f= divide (f, c_f, as);
-
-  f *= c_gcd;
-  f /= vcontent (f, Variable (v+1));
-
-  return f;
-}
-
+
+
Index: factory/libfac/charset/csutil.h
===================================================================
--- factory/libfac/charset/csutil.h	(revision 24819474a270d2480ed2d9f37a85ab3f5eafee1c)
+++ factory/libfac/charset/csutil.h	(revision 5355082c234c56e9e117c19600ae46cf6324f5ca)
@@ -52,7 +52,5 @@
 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 );
Index: factory/libfac/factor/helpstuff.cc
===================================================================
--- factory/libfac/factor/helpstuff.cc	(revision 24819474a270d2480ed2d9f37a85ab3f5eafee1c)
+++ factory/libfac/factor/helpstuff.cc	(revision 5355082c234c56e9e117c19600ae46cf6324f5ca)
@@ -21,44 +21,2 @@
 }
 
-// Now some procedures used in SqrFree and in Factor
-///////////////////////////////////////////////////////////////
-///////////////////////////////////////////////////////////////
-// We have to include a version of <CFFList>.append(CFFactor)//
-// and Union( CFFList, CFFList)                              //
-// because we have to look for multiplicities in SqrFree.    //
-// e.g.: SqrFree( f^3 ) with char <> 3                       //
-///////////////////////////////////////////////////////////////
-CFFList
-myappend( const CFFList & Inputlist, const CFFactor & TheFactor)
-{
-  CFFList Outputlist ;
-  CFFactor copy;
-  CFFListIterator i;
-  int exp=0;
-
-  for ( i=Inputlist ; i.hasItem() ; i++ )
-  {
-    copy = i.getItem();
-    if ( copy.factor() == TheFactor.factor() )
-      exp += copy.exp();
-    else
-      Outputlist.append(copy);
-  }
-  Outputlist.append( CFFactor(TheFactor.factor(), exp + TheFactor.exp()));
-  return Outputlist;
-}
-
-CFFList
-myUnion(const CFFList & Inputlist1,const CFFList & Inputlist2)
-{
-  CFFList Outputlist;
-  CFFListIterator i;
-
-  for ( i=Inputlist1 ; i.hasItem() ; i++ )
-    Outputlist = myappend(Outputlist, i.getItem() );
-  for ( i=Inputlist2 ; i.hasItem() ; i++ )
-    Outputlist = myappend(Outputlist, i.getItem() );
-
-  return Outputlist;
-}
-
Index: factory/libfac/factor/helpstuff.h
===================================================================
--- factory/libfac/factor/helpstuff.h	(revision 24819474a270d2480ed2d9f37a85ab3f5eafee1c)
+++ factory/libfac/factor/helpstuff.h	(revision 5355082c234c56e9e117c19600ae46cf6324f5ca)
@@ -9,6 +9,4 @@
 // Now some procedures used in SqrFree and in Factor
 ///////////////////////////////////////////////////////////////
-CFFList myappend( const CFFList & Inputlist, const CFFactor & TheFactor) ;
-CFFList myUnion(const CFFList & Inputlist1,const CFFList & Inputlist2);
 inline int min ( const int a, const int b ){
   return (a<=b ? a:b);
