cfEzgcd.h File Reference

Extended Zassenhaus GCD over finite fields and Z. More...

Go to the source code of this file.

Functions
CanonicalForm	EZGCD_P (const CanonicalForm &A, const CanonicalForm &B)
	Extended Zassenhaus GCD for finite fields. In case things become too dense we switch to a modular algorithm. More...
CanonicalForm	ezgcd (const CanonicalForm &f, const CanonicalForm &g)
	Extended Zassenhaus GCD over Z. In case things become too dense we switch to a modular algorithm. More...

Detailed Description

Extended Zassenhaus GCD over finite fields and Z.

Definition in file cfEzgcd.h.

Function Documentation

ezgcd()

CanonicalForm ezgcd	(	const CanonicalForm &	f,
		const CanonicalForm &	g
	)

Extended Zassenhaus GCD over Z. In case things become too dense we switch to a modular algorithm.

Definition at line 861 of file cfEzgcd.cc.

{
  REvaluation b;
  return ezgcd( FF, GG, b, false );
}

EZGCD_P()

CanonicalForm EZGCD_P	(	const CanonicalForm &	A,
		const CanonicalForm &	B
	)

Extended Zassenhaus GCD for finite fields. In case things become too dense we switch to a modular algorithm.

Definition at line 876 of file cfEzgcd.cc.

{
  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
  else if (FF.isZero() && GG.isZero()) return FF.genOne();
  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
  if (FF == GG) return FF/Lc(FF);
 
  int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
  Variable a, oldA;
  int sizeF= size (FF);
  int sizeG= size (GG);
 
  if (sizeF==1)
  {
    return gcd_mon( FF, GG );
  }
  else if (sizeG==1)
  {
    return gcd_mon( GG, FF );
  }
 
  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
  {
    if (hasFirstAlgVar (FF, a) || hasFirstAlgVar (GG, a))
      return modGCDFq (FF, GG, a);
    else if (CFFactory::gettype() == GaloisFieldDomain)
      return modGCDGF (FF, GG);
    else
      return modGCDFp (FF, GG);
  }
 
  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
                lcD;
  CFArray DD( 1, 2 ), lcDD( 1, 2 );
  int degF, degG, delta, count;
  int maxeval;
  int ch=getCharacteristic();
  maxeval= tmin((ch/
                (int)(ilog2(ch)*log2exp))*2, maxNumEval);
  count= 0; // number of eval. used
  REvaluation b, bt;
  int gcdfound = 0;
  Variable x = Variable(1);
 
  F= FF;
  G= GG;
 
  CFMap M,N;
  int smallestDegLev;
  TIMING_DEFINE(ez_p_compress);
  TIMING_START (ez_p_compress);
  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
 
  if (best_level == 0) return G.genOne();
 
  F= M (F);
  G= M (G);
  TIMING_END_AND_PRINT (ez_p_compress, "time for compression in EZ_P: ")
 
  TIMING_DEFINE (ez_p_content)
  TIMING_START (ez_p_content)
  f = content( F, x ); g = content( G, x );
  d = gcd( f, g );
  F /= f; G /= g;
  TIMING_END_AND_PRINT (ez_p_content, "time to extract content in EZ_P: ")
 
  if( F.isUnivariate() && G.isUnivariate() )
  {
    if( F.mvar() == G.mvar() )
      d *= gcd( F, G );
    else
      return N (d);
    return N (d);
  }
  if ( F.isUnivariate())
  {
    g= content (G,G.mvar());
    return N(d*gcd(F,g));
  }
  if ( G.isUnivariate())
  {
    f= content (F,F.mvar());
    return N(d*gcd(G,f));
  }
 
  maxNumVars= tmax (getNumVars (F), getNumVars (G));
  sizeF= size (F);
  sizeG= size (G);
 
  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
  {
    if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
      return N (d*modGCDFq (F, G, a));
    else if (CFFactory::gettype() == GaloisFieldDomain)
      return N (d*modGCDGF (F, G));
    else
      return N (d*modGCDFp (F, G));
  }
 
  int dummy= 0;
  if( gcd_test_one( F, G, false, dummy ) )
  {
    return N (d);
  }
 
  bool passToGF= false;
  bool extOfExt= false;
  int p= getCharacteristic();
  bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
  int k= 1;
  CanonicalForm primElem, imPrimElem;
  CFList source, dest;
  if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
  {
    if (p == 2)
      setCharacteristic (2, 12, 'Z');
    else if (p == 3)
      setCharacteristic (3, 4, 'Z');
    else if (p == 5 || p == 7)
      setCharacteristic (p, 3, 'Z');
    else
      setCharacteristic (p, 2, 'Z');
    passToGF= true;
    F= F.mapinto();
    G= G.mapinto();
    maxeval= 2*ipower (p, getGFDegree());
  }
  else if (CFFactory::gettype() == GaloisFieldDomain &&
           ipower (p , getGFDegree()) < 50)
  {
    k= getGFDegree();
    if (ipower (p, 2*k) > 50)
      setCharacteristic (p, 2*k, gf_name);
    else
      setCharacteristic (p, 3*k, gf_name);
    F= GFMapUp (F, k);
    G= GFMapUp (G, k);
    maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
  }
  else if (p < 50 && algExtension && CFFactory::gettype() != GaloisFieldDomain)
  {
    int d= degree (getMipo (a));
    oldA= a;
    Variable v2;
    if (p == 2 && d < 6)
    {
      bool primFail= false;
      Variable vBuf;
      primElem= primitiveElement (a, vBuf, primFail);
      ASSERT (!primFail, "failure in integer factorizer");
      if (d < 3)
      {
        #ifdef HAVE_FLINT
        nmod_poly_t Irredpoly;
        nmod_poly_init(Irredpoly,p);
        nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, 3*d+1);
        CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
        nmod_poly_clear(Irredpoly);
        #elif defined(HAVE_NTL)
        if (fac_NTL_char != p)
        {
          fac_NTL_char= p;
          zz_p::init (p);
        }
        zz_pX NTLIrredpoly;
        BuildIrred (NTLIrredpoly, d*3);
        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
        #else
        factoryError("NTL/FLINT missing: EZGCD_P");
        #endif
        v2= rootOf (newMipo);
      }
      else
      {
        #ifdef HAVE_FLINT
        nmod_poly_t Irredpoly;
        nmod_poly_init(Irredpoly,p);
        nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, 2*d+1);
        CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
        nmod_poly_clear(Irredpoly);
        #elif defined(HAVE_NTL)
        if (fac_NTL_char != p)
        {
          fac_NTL_char= p;
          zz_p::init (p);
        }
        zz_pX NTLIrredpoly;
        BuildIrred (NTLIrredpoly, d*2);
        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
        #else
        factoryError("NTL/FLINT missing: EZGCD_P");
        #endif
        v2= rootOf (newMipo);
      }
      imPrimElem= mapPrimElem (primElem, a, v2);
      extOfExt= true;
    }
    else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
    {
      bool primFail= false;
      Variable vBuf;
      primElem= primitiveElement (a, vBuf, primFail);
      ASSERT (!primFail, "failure in integer factorizer");
      #ifdef HAVE_FLINT
      nmod_poly_t Irredpoly;
      nmod_poly_init(Irredpoly,p);
      nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, 2*d+1);
      CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
      nmod_poly_clear(Irredpoly);
      #elif defined(HAVE_NTL)
      if (fac_NTL_char != p)
      {
        fac_NTL_char= p;
        zz_p::init (p);
      }
      zz_pX NTLIrredpoly;
      BuildIrred (NTLIrredpoly, d*2);
      CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
      #else
      factoryError("NTL/FLINT missing: EZGCD_P");
      #endif
      v2= rootOf (newMipo);
      imPrimElem= mapPrimElem (primElem, a, v2);
      extOfExt= true;
    }
    if (extOfExt)
    {
      maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
      F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
      G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
      a= v2;
    }
  }
 
  lcF = LC( F, x ); lcG = LC( G, x );
  lcD = gcd( lcF, lcG );
 
  delta = 0;
  degF = degree( F, x ); degG = degree( G, x );
 
  if (algExtension)
    b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
  else
  { // both not in extension given by algebraic variable
    if (CFFactory::gettype() != GaloisFieldDomain)
      b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
    else
      b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
  }
 
  CanonicalForm cand, contcand;
  CanonicalForm result;
  int o, t;
  o= 0;
  t= 1;
  int goodPointCount= 0;
  TIMING_DEFINE(ez_p_eval);
  while( !gcdfound )
  {
    TIMING_START (ez_p_eval);
    if( !findeval( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
         maxeval/maxNumVars, t ))
    { // too many eval. used --> try another method
      Off (SW_USE_EZGCD_P);
      result= gcd (F,G);
      On (SW_USE_EZGCD_P);
      if (passToGF)
      {
        CanonicalForm mipo= gf_mipo;
        setCharacteristic (p);
        Variable alpha= rootOf (mipo.mapinto());
        result= GF2FalphaRep (result, alpha);
        prune (alpha);
      }
      if (k > 1)
      {
        result= GFMapDown (result, k);
        setCharacteristic (p, k, gf_name);
      }
      if (extOfExt)
      {
        result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
        prune1 (oldA);
      }
      return N (d*result);
    }
    TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P1: ");
    delta = degree( Db );
    if (delta == degF)
    {
      if (degF <= degG && fdivides (F, G))
      {
        if (passToGF)
        {
          CanonicalForm mipo= gf_mipo;
          setCharacteristic (p);
          Variable alpha= rootOf (mipo.mapinto());
          F= GF2FalphaRep (F, alpha);
          prune (alpha);
        }
        if (k > 1)
        {
          F= GFMapDown (F, k);
          setCharacteristic (p, k, gf_name);
        }
        if (extOfExt)
        {
          F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
          prune1 (oldA);
        }
        return N (d*F);
      }
      else
        delta--;
    }
    else if (delta == degG)
    {
      if (degG <= degF && fdivides (G, F))
      {
        if (passToGF)
        {
          CanonicalForm mipo= gf_mipo;
          setCharacteristic (p);
          Variable alpha= rootOf (mipo.mapinto());
          G= GF2FalphaRep (G, alpha);
          prune (alpha);
        }
        if (k > 1)
        {
          G= GFMapDown (G, k);
          setCharacteristic (p, k, gf_name);
        }
        if (extOfExt)
        {
          G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
          prune1 (oldA);
        }
        return N (d*G);
      }
      else
        delta--;
    }
    if( delta == 0 )
    {
      if (passToGF)
        setCharacteristic (p);
      if (k > 1)
        setCharacteristic (p, k, gf_name);
      return N (d);
    }
    while( true )
    {
      bt = b;
      TIMING_START (ez_p_eval);
      if( !findeval(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
           maxeval/maxNumVars, t ))
      { // too many eval. used --> try another method
        Off (SW_USE_EZGCD_P);
        result= gcd (F,G);
        On (SW_USE_EZGCD_P);
        if (passToGF)
        {
          CanonicalForm mipo= gf_mipo;
          setCharacteristic (p);
          Variable alpha= rootOf (mipo.mapinto());
          result= GF2FalphaRep (result, alpha);
          prune (alpha);
        }
        if (k > 1)
        {
          result= GFMapDown (result, k);
          setCharacteristic (p, k, gf_name);
        }
        if (extOfExt)
        {
          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
          prune1 (oldA);
        }
        return N (d*result);
      }
      TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P2: ");
      int dd = degree( Dbt );
      if( dd == 0 )
      {
        if (passToGF)
          setCharacteristic (p);
        if (k > 1)
          setCharacteristic (p, k, gf_name);
        return N (d);
      }
      if( dd == delta )
      {
        goodPointCount++;
        if (goodPointCount == 5)
          break;
      }
      if( dd < delta )
      {
        goodPointCount= 0;
        delta = dd;
        b = bt;
        Db = Dbt; Fb = Fbt; Gb = Gbt;
      }
      if (delta == degF)
      {
        if (degF <= degG && fdivides (F, G))
        {
          if (passToGF)
          {
            CanonicalForm mipo= gf_mipo;
            setCharacteristic (p);
            Variable alpha= rootOf (mipo.mapinto());
            F= GF2FalphaRep (F, alpha);
            prune (alpha);
          }
          if (k > 1)
          {
            F= GFMapDown (F, k);
            setCharacteristic (p, k, gf_name);
          }
          if (extOfExt)
          {
            F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
            prune1 (oldA);
          }
          return N (d*F);
        }
        else
          delta--;
      }
      else if (delta == degG)
      {
        if (degG <= degF && fdivides (G, F))
        {
          if (passToGF)
          {
            CanonicalForm mipo= gf_mipo;
            setCharacteristic (p);
            Variable alpha= rootOf (mipo.mapinto());
            G= GF2FalphaRep (G, alpha);
            prune (alpha);
          }
          if (k > 1)
          {
            G= GFMapDown (G, k);
            setCharacteristic (p, k, gf_name);
          }
          if (extOfExt)
          {
            G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
            prune1 (oldA);
          }
          return N (d*G);
        }
        else
          delta--;
      }
      if( delta == 0 )
      {
        if (passToGF)
          setCharacteristic (p);
        if (k > 1)
          setCharacteristic (p, k, gf_name);
        return N (d);
      }
    }
    if( delta != degF && delta != degG )
    {
      bool B_is_F;
      CanonicalForm xxx1, xxx2;
      CanonicalForm buf;
      DD[1] = Fb / Db;
      buf= Gb/Db;
      xxx1 = gcd( DD[1], Db );
      xxx2 = gcd( buf, Db );
      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
          (size (F) <= size (G)))
          || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
      {
        B = F;
        DD[2] = Db;
        lcDD[1] = lcF;
        lcDD[2] = lcD;
        B_is_F = true;
      }
      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
               (size (G) < size (F)))
               || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
      {
        DD[1] = buf;
        B = G;
        DD[2] = Db;
        lcDD[1] = lcG;
        lcDD[2] = lcD;
        B_is_F = false;
      }
      else // special case handling
      {
        Off (SW_USE_EZGCD_P);
        result= gcd (F,G);
        On (SW_USE_EZGCD_P);
        if (passToGF)
        {
          CanonicalForm mipo= gf_mipo;
          setCharacteristic (p);
          Variable alpha= rootOf (mipo.mapinto());
          result= GF2FalphaRep (result, alpha);
          prune (alpha);
        }
        if (k > 1)
        {
          result= GFMapDown (result, k);
          setCharacteristic (p, k, gf_name);
        }
        if (extOfExt)
        {
          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
          prune1 (oldA);
        }
        return N (d*result);
      }
      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
 
      if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
      {
        if (algExtension)
        {
          result= modGCDFq (F, G, a);
          if (extOfExt)
          {
            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
            prune1 (oldA);
          }
          return N (d*result);
        }
        if (CFFactory::gettype() == GaloisFieldDomain)
        {
          result= modGCDGF (F, G);
          if (passToGF)
          {
            CanonicalForm mipo= gf_mipo;
            setCharacteristic (p);
            Variable alpha= rootOf (mipo.mapinto());
            result= GF2FalphaRep (result, alpha);
            prune (alpha);
          }
          if (k > 1)
          {
            result= GFMapDown (result, k);
            setCharacteristic (p, k, gf_name);
          }
          return N (d*result);
        }
        else
          return N (d*modGCDFp (F,G));
      }
 
      TIMING_DEFINE(ez_p_hensel_lift);
      TIMING_START (ez_p_hensel_lift);
      gcdfound= Hensel (B*lcD, DD, b, lcDD);
      TIMING_END_AND_PRINT (ez_p_hensel_lift, "time for Hensel lift in EZ_P: ");
 
      if (gcdfound == -1) //things became dense
      {
        if (algExtension)
        {
          result= modGCDFq (F, G, a);
          if (extOfExt)
          {
            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
            prune1 (oldA);
          }
          return N (d*result);
        }
        if (CFFactory::gettype() == GaloisFieldDomain)
        {
          result= modGCDGF (F, G);
          if (passToGF)
          {
            CanonicalForm mipo= gf_mipo;
            setCharacteristic (p);
            Variable alpha= rootOf (mipo.mapinto());
            result= GF2FalphaRep (result, alpha);
            prune (alpha);
          }
          if (k > 1)
          {
            result= GFMapDown (result, k);
            setCharacteristic (p, k, gf_name);
          }
          return N (d*result);
        }
        else
        {
          if (p >= cf_getBigPrime(0))
            return N (d*sparseGCDFp (F,G));
          else
            return N (d*modGCDFp (F,G));
        }
      }
 
      if (gcdfound == 1)
      {
        TIMING_DEFINE(termination_test);
        TIMING_START (termination_test);
        contcand= content (DD[2], Variable (1));
        cand = DD[2] / contcand;
        if (B_is_F)
          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
        else
          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
        TIMING_END_AND_PRINT (termination_test,
                              "time for termination test EZ_P: ");
 
        if (passToGF && gcdfound)
        {
          CanonicalForm mipo= gf_mipo;
          setCharacteristic (p);
          Variable alpha= rootOf (mipo.mapinto());
          cand= GF2FalphaRep (cand, alpha);
          prune (alpha);
        }
        if (k > 1 && gcdfound)
        {
          cand= GFMapDown (cand, k);
          setCharacteristic (p, k, gf_name);
        }
        if (extOfExt && gcdfound)
        {
          cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
          prune1 (oldA);
        }
      }
    }
    delta--;
    goodPointCount= 0;
  }
  return N (d*cand);
}