cfModResultant.cc File Reference

modular resultant algorithm More...

#include "config.h"
#include "cf_assert.h"
#include "timing.h"
#include "cfModResultant.h"
#include "cf_primes.h"
#include "templates/ftmpl_functions.h"
#include "cf_map.h"
#include "cf_algorithm.h"
#include "cf_iter.h"
#include "cf_irred.h"
#include "cf_generator.h"
#include "cf_random.h"
#include "cf_map_ext.h"
#include "facFqBivarUtil.h"
#include "NTLconvert.h"
#include "FLINTconvert.h"

Go to the source code of this file.

Functions
	TIMING_DEFINE_PRINT (fac_resultant_p) static inline void myCompress(const CanonicalForm &F
	for (int i=0;i<=n;i++) degsf[i]
	if (x.level() !=1)
	DELETE_ARRAY (degsg)
static CanonicalForm	oneNorm (const CanonicalForm &F)
static CanonicalForm	uniResultant (const CanonicalForm &F, const CanonicalForm &G)
static void	evalPoint (const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &FEval, CanonicalForm &GEval, CFGenerator &evalPoint)
static CanonicalForm	newtonInterp (const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x)
CanonicalForm	resultantFp (const CanonicalForm &A, const CanonicalForm &B, const Variable &x, bool prob)
	modular resultant algorihtm over Fp More...
static CanonicalForm	balanceUni (const CanonicalForm &f, const CanonicalForm &q)
static CanonicalForm	symmetricRemainder (const CanonicalForm &f, const CanonicalForm &q)
CanonicalForm	resultantZ (const CanonicalForm &A, const CanonicalForm &B, const Variable &x, bool prob)
	modular resultant algorihtm over Z More...

Variables
const CanonicalForm &	G
const CanonicalForm CFMap &	M
const CanonicalForm CFMap CFMap &	N
const CanonicalForm CFMap CFMap const Variable &	x
int *	degsf = NEW_ARRAY(int,n + 1)
int *	degsg = NEW_ARRAY(int,n + 1)
int	both_non_zero = 0
int	f_zero = 0
int	g_zero = 0
int	both_zero = 0
int	degsfx
int	degsgx
int	Flevel =F.level()
int	Glevel =G.level()
	else

Detailed Description

modular resultant algorithm

Author

Martin Lee

Definition in file cfModResultant.cc.

Function Documentation

balanceUni()

static CanonicalForm balanceUni ( const CanonicalForm &  f, const CanonicalForm &  q  )

inlinestatic

Definition at line 526 of file cfModResultant.cc.

{
    Variable x = f.mvar();
    CanonicalForm result = 0, qh = q / 2;
    CanonicalForm c;
    CFIterator i;
    for ( i = f; i.hasTerms(); i++ ) {
        c = mod( i.coeff(), q );
        if ( c > qh )
            result += power( x, i.exp() ) * (c - q);
        else
            result += power( x, i.exp() ) * c;
    }
    return result;
}

DELETE_ARRAY()

DELETE_ARRAY

(

degsg

)

evalPoint()

static void evalPoint ( const CanonicalForm &  F, const CanonicalForm &  G, CanonicalForm &  FEval, CanonicalForm &  GEval, CFGenerator &  evalPoint  )

inlinestatic

Definition at line 310 of file cfModResultant.cc.

{
  int degF, degG;
  Variable x= Variable (1);
  degF= degree (F, x);
  degG= degree (G, x);
  do
  {
    if (!evalPoint.hasItems())
      break;
    FEval= F (evalPoint.item(), 2);
    GEval= G (evalPoint.item(), 2);
    if (degree (FEval, 1) < degF || degree (GEval, 1) < degG)
    {
      evalPoint.next();
      continue;
    }
    else
      return;
  }
  while (evalPoint.hasItems());
}

for()

for

(

int

i = 0;i<=n;i++

)

Definition at line 72 of file cfEzgcd.cc.

  {
    if (degsf[i] != 0 && degsg[i] != 0)
    {
      both_non_zero++;
      continue;
    }
    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
    {
      f_zero++;
      continue;
    }
    if (degsg[i] == 0 && degsf[i] && i <= F.level())
    {
      g_zero++;
      continue;
    }
  }

if()

if

(

x.level() !

= 1

)

Definition at line 66 of file cfModResultant.cc.

  {
    int xlevel= x.level();
 
    for (int i= 1; i <= n; i++)
    {
      if (degsf[i] != 0 && degsg[i] != 0)
      {
        both_non_zero++;
        continue;
      }
      if (degsf[i] == 0 && degsg[i] != 0 && i <= Glevel)
      {
        f_zero++;
        continue;
      }
      if (degsg[i] == 0 && degsf[i] && i <= Flevel)
      {
        g_zero++;
        continue;
      }
    }
 
    M.newpair (Variable (xlevel), Variable (1));
    N.newpair (Variable (1), Variable (xlevel));
    degsfx= degsf [xlevel];
    degsgx= degsg [xlevel];
    degsf [xlevel]= 0;
    degsg [xlevel]= 0;
    if ((getNumVars (F) == 2 && getNumVars (G) == 1) ||
        (getNumVars (G) == 2 && getNumVars (F) == 1) ||
        (getNumVars (F) == 2 && getNumVars (F) == getNumVars (G)
         && getVars (F) == getVars (G)))
    {
      int pos= 2;
      for (int i= 1; i <= n; i++)
      {
        if (i != xlevel)
        {
          if (i != pos && (degsf[i] != 0 || degsg [i] != 0))
          {
            M.newpair (Variable (i), Variable (pos));
            N.newpair (Variable (pos), Variable (i));
            pos++;
          }
        }
      }
      DELETE_ARRAY(degsf);
      DELETE_ARRAY(degsg);
      return;
    }
 
    if (both_non_zero == 0)
    {
      DELETE_ARRAY(degsf);
      DELETE_ARRAY(degsg);
      return;
    }
 
    // map Variables which do not occur in both polynomials to higher levels
    int k= 1;
    int l= 1;
    for (int i= 1; i <= n; i++)
    {
      if (i == xlevel)
        continue;
      if (degsf[i] != 0 && degsg[i] == 0 && i <= Flevel)
      {
        if (k + both_non_zero != i)
        {
          M.newpair (Variable (i), Variable (k + both_non_zero));
          N.newpair (Variable (k + both_non_zero), Variable (i));
        }
        k++;
      }
      if (degsf[i] == 0 && degsg[i] != 0 && i <= Glevel)
      {
        if (l + g_zero + both_non_zero != i)
        {
          M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
          N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
        }
        l++;
      }
    }
 
    int m= n;
    int min_max_deg;
    k= both_non_zero;
    l= 0;
    int i= 1;
    while (k > 0)
    {
      if (degsf [i] != 0 && degsg [i] != 0)
        min_max_deg= degsgx*degsf[i] + degsfx*degsg[i];
      else
        min_max_deg= 0;
      while (min_max_deg == 0 && i < m)
      {
        i++;
        if (degsf [i] != 0 && degsg [i] != 0)
          min_max_deg= degsgx*degsf[i] + degsfx*degsg[i];
        else
          min_max_deg= 0;
      }
      for (int j= i + 1; j <= m; j++)
      {
        if (degsgx*degsf[j] + degsfx*degsg[j] <= min_max_deg &&
            degsf[j] != 0 && degsg [j] != 0)
        {
          min_max_deg= degsgx*degsf[j] + degsfx*degsg[j];
          l= j;
        }
      }
      if (l != 0)
      {
        if (l != k && l != xlevel && k != 1)
        {
          M.newpair (Variable (l), Variable(k));
          N.newpair (Variable (k), Variable(l));
          degsf[l]= 0;
          degsg[l]= 0;
          l= 0;
        }
        else if (l < m + 1)
        {
          degsf[l]= 0;
          degsg[l]= 0;
          l= 0;
        }
      }
      else if (l == 0)
      {
        if (i != k && i != xlevel && k != 1)
        {
          M.newpair (Variable (i), Variable (k));
          N.newpair (Variable (k), Variable (i));
          degsf[i]= 0;
          degsg[i]= 0;
        }
        else if (i < m + 1)
        {
          degsf[i]= 0;
          degsg[i]= 0;
        }
        i++;
      }
      k--;
    }
  }

newtonInterp()

static CanonicalForm newtonInterp ( const CanonicalForm &  alpha, const CanonicalForm &  u, const CanonicalForm &  newtonPoly, const CanonicalForm &  oldInterPoly, const Variable &  x  )

inlinestatic

Definition at line 336 of file cfModResultant.cc.

{
  CanonicalForm interPoly;
 
  interPoly= oldInterPoly+((u - oldInterPoly (alpha, x))/newtonPoly (alpha, x))
                            *newtonPoly;
  return interPoly;
}

oneNorm()

static CanonicalForm oneNorm ( const CanonicalForm &  F )

inlinestatic

Definition at line 242 of file cfModResultant.cc.

{
  if (F.inZ())
    return abs (F);
 
  CanonicalForm result= 0;
  for (CFIterator i= F; i.hasTerms(); i++)
    result += oneNorm (i.coeff());
  return result;
}

resultantFp()

CanonicalForm resultantFp	(	const CanonicalForm &	A,
		const CanonicalForm &	B,
		const Variable &	x,
		bool	prob = true
	)

modular resultant algorihtm over Fp

Returns

resultantFp returns the resultant of A and B wrt. x

Parameters

[in]	A	some poly
[in]	B	some poly
[in]	x	some polynomial variable
[in]	prob	if true use probabilistic algorithm

Definition at line 348 of file cfModResultant.cc.

{
  ASSERT (getCharacteristic() > 0, "characteristic > 0 expected");
 
  if (A.isZero() || B.isZero())
    return 0;
 
  int degAx= degree (A, x);
  int degBx= degree (B, x);
  if (A.level() < x.level())
    return power (A, degBx);
  if (B.level() < x.level())
    return power (B, degAx);
 
  if (degAx == 0)
    return power (A, degBx);
  else if (degBx == 0)
    return power (B, degAx);
 
  if (A.isUnivariate() && B.isUnivariate() && A.level() == B.level())
    return uniResultant (A, B);
 
  CanonicalForm F= A;
  CanonicalForm G= B;
 
  CFMap M, N;
  myCompress (F, G, M, N, x);
 
  F= M (F);
  G= M (G);
 
  Variable y= Variable (2);
 
  CanonicalForm GEval, FEval, recResult, H;
  CanonicalForm newtonPoly= 1;
  CanonicalForm modResult= 0;
 
  Variable z= Variable (1);
  int bound= degAx*degree (G, 2) + degree (F, 2)*degBx;
 
  int p= getCharacteristic();
  CanonicalForm minpoly;
  Variable alpha= Variable (tmax (F.level(), G.level()) + 1);
  bool algExt= hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha);
  CFGenerator * gen;
  bool extOfExt= false;
  Variable v= alpha;
  CanonicalForm primElemAlpha, imPrimElemAlpha;
  CFList source,dest;
  if (!algExt && (p < (1 << 28)))
  {
    // pass to an extension of size at least 2^29
    // for very very large input that is maybe too small though
    int deg= ceil (29.0*((double) log (2)/log (p)))+1;
    minpoly= randomIrredpoly (deg, z);
    alpha= rootOf (minpoly);
    AlgExtGenerator AlgExtGen (alpha);
    gen= AlgExtGen.clone();
    for (int i= 0; i < p; i++) // skip values from the prime field
      (*gen).next();
  }
  else if (!algExt)
  {
    FFGenerator FFGen;
    gen= FFGen.clone();
  }
  else
  {
    int deg= ceil (29.0*((double) log (2)/log (p)));
    if (degree (getMipo (alpha)) < deg)
    {
      mpz_t field_size;
      mpz_init (field_size);
      mpz_ui_pow_ui (field_size, p,
                 deg + degree (getMipo (alpha)) - deg%degree (getMipo (alpha)));
 
      // field_size needs to fit in an int because of mapUp, mapDown, length of lists etc.
      if (mpz_fits_sint_p (field_size))
      {
        minpoly= randomIrredpoly (deg + degree (getMipo (alpha))
                                  - deg%degree (getMipo (alpha)), z);
        v= rootOf (minpoly);
        Variable V_buf2;
        bool primFail= false;
        extOfExt= true;
        primElemAlpha= primitiveElement (alpha, V_buf2, primFail);
        ASSERT (!primFail, "failure in integer factorizer");
        if (primFail)
          ; //ERROR
        else
          imPrimElemAlpha= mapPrimElem (primElemAlpha, alpha, v);
        F= mapUp (F, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
        G= mapUp (G, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
      }
      else
      {
        deg= deg - deg % degree (getMipo (alpha));
        mpz_ui_pow_ui (field_size, p, deg);
        while (deg / degree (getMipo (alpha)) >= 2 && !mpz_fits_sint_p (field_size))
        {
          deg -= degree (getMipo (alpha));
          mpz_ui_pow_ui (field_size, p, deg);
        }
        if (deg != degree (getMipo (alpha)))
        {
           minpoly= randomIrredpoly (deg, z);
           v= rootOf (minpoly);
           Variable V_buf2;
           bool primFail= false;
           extOfExt= true;
           primElemAlpha= primitiveElement (alpha, V_buf2, primFail);
           ASSERT (!primFail, "failure in integer factorizer");
           if (primFail)
             ; //ERROR
           else
             imPrimElemAlpha= mapPrimElem (primElemAlpha, alpha, v);
           F= mapUp (F, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
           G= mapUp (G, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
        }
      }
      mpz_clear (field_size);
    }
    AlgExtGenerator AlgExtGen (v);
    gen= AlgExtGen.clone();
    for (int i= 0; i < p; i++)
      (*gen).next();
  }
  int count= 0;
  int equalCount= 0;
  CanonicalForm point;
  do
  {
    evalPoint (F, G, FEval, GEval, *gen);
 
    recResult= resultantFp (FEval, GEval, z, prob);
 
    H= newtonInterp ((*gen).item(), recResult, newtonPoly, modResult, y);
 
    if (H == modResult)
      equalCount++;
    else
      equalCount= 0;
 
    count++;
    if (count > bound || (prob && equalCount == 2 && !H.inCoeffDomain()))
    {
      if (!algExt && degree (H, alpha) <= 0)
        break;
      else if (algExt)
      {
        if (extOfExt && !isInExtension (H, imPrimElemAlpha, 1, primElemAlpha,
                                        dest, source))
        {
          H= mapDown (H, primElemAlpha, imPrimElemAlpha, alpha, dest, source);
          prune (v);
          break;
        }
        else if (!extOfExt)
          break;
      }
    }
 
    modResult= H;
    newtonPoly *= (y - (*gen).item());
    if ((*gen).hasItems())
        (*gen).next();
    else
      STICKYASSERT (0, "out of evaluation points");
  } while (1);
 
  delete gen;
 
  return N (H);
}

resultantZ()

CanonicalForm resultantZ	(	const CanonicalForm &	A,
		const CanonicalForm &	B,
		const Variable &	x,
		bool	prob = true
	)

modular resultant algorihtm over Z

Returns

resultantZ returns the resultant of A and B wrt. x

Parameters

[in]	A	some poly
[in]	B	some poly
[in]	x	some polynomial variable
[in]	prob	if true use probabilistic algorithm

Definition at line 559 of file cfModResultant.cc.

{
  ASSERT (getCharacteristic() == 0, "characteristic > 0 expected");
#ifndef NOASSERT
  bool isRat= isOn (SW_RATIONAL);
  On (SW_RATIONAL);
  ASSERT (bCommonDen (A).isOne(), "input A is rational");
  ASSERT (bCommonDen (B).isOne(), "input B is rational");
  if (!isRat)
    Off (SW_RATIONAL);
#endif
 
  int degAx= degree (A, x);
  int degBx= degree (B, x);
  if (A.level() < x.level())
    return power (A, degBx);
  if (B.level() < x.level())
    return power (B, degAx);
 
  if (degAx == 0)
    return power (A, degBx);
  else if (degBx == 0)
    return power (B, degAx);
 
  CanonicalForm F= A;
  CanonicalForm G= B;
 
  Variable X= x;
  if (F.level() != x.level() || G.level() != x.level())
  {
    if (F.level() > G.level())
      X= F.mvar();
    else
      X= G.mvar();
    F= swapvar (F, X, x);
    G= swapvar (G, X, x);
  }
 
  // now X is the main variable
 
  CanonicalForm d= 0;
  CanonicalForm dd= 0;
  CanonicalForm buf;
  for (CFIterator i= F; i.hasTerms(); i++)
  {
    buf= oneNorm (i.coeff());
    d= (buf > d) ? buf : d;
  }
  CanonicalForm e= 0, ee= 0;
  for (CFIterator i= G; i.hasTerms(); i++)
  {
    buf= oneNorm (i.coeff());
    e= (buf > e) ? buf : e;
  }
  d= power (d, degBx);
  e= power (e, degAx);
  CanonicalForm bound= 1;
  for (int i= degBx + degAx; i > 1; i--)
    bound *= i;
  bound *= d*e;
  bound *= 2;
 
  bool onRational= isOn (SW_RATIONAL);
  if (onRational)
    Off (SW_RATIONAL);
  int i = cf_getNumBigPrimes() - 1;
  int p;
  CanonicalForm l= lc (F)*lc(G);
  CanonicalForm resultModP, q (0), newResult, newQ;
  CanonicalForm result;
  int equalCount= 0;
  CanonicalForm test, newTest;
  int count= 0;
  do
  {
    p = cf_getBigPrime( i );
    i--;
    while ( i >= 0 && mod( l, p ) == 0)
    {
      p = cf_getBigPrime( i );
      i--;
    }
 
    if (i <= 0)
      return resultant (A, B, x);
 
    setCharacteristic (p);
 
    TIMING_START (fac_resultant_p);
    resultModP= resultantFp (mapinto (F), mapinto (G), X, prob);
    TIMING_END_AND_PRINT (fac_resultant_p, "time to compute resultant mod p: ");
 
    setCharacteristic (0);
 
    count++;
    if ( q.isZero() )
    {
      result= mapinto(resultModP);
      q= p;
    }
    else
    {
      chineseRemainder( result, q, mapinto (resultModP), p, newResult, newQ );
      q= newQ;
      result= newResult;
      test= symmetricRemainder (result,q);
      if (test != newTest)
      {
        newTest= test;
        equalCount= 0;
      }
      else
        equalCount++;
      if (newQ > bound || (prob && equalCount == 2))
      {
        result= test;
        break;
      }
    }
  } while (1);
 
  if (onRational)
    On (SW_RATIONAL);
  return swapvar (result, X, x);
}

symmetricRemainder()

static CanonicalForm symmetricRemainder ( const CanonicalForm &  f, const CanonicalForm &  q  )

inlinestatic

Definition at line 544 of file cfModResultant.cc.

{
  CanonicalForm result= 0;
  if (f.isUnivariate() || f.inCoeffDomain())
    return balanceUni (f, q);
  else
  {
    Variable x= f.mvar();
    for (CFIterator i= f; i.hasTerms(); i++)
      result += power (x, i.exp())*symmetricRemainder (i.coeff(), q);
  }
  return result;
}

TIMING_DEFINE_PRINT()

TIMING_DEFINE_PRINT

(

fac_resultant_p

)

const &

uniResultant()

static CanonicalForm uniResultant ( const CanonicalForm &  F, const CanonicalForm &  G  )

inlinestatic

Definition at line 255 of file cfModResultant.cc.

{
  ASSERT (getCharacteristic() > 0, "characteristic > 0 expected");
  if (F.inCoeffDomain() && G.inCoeffDomain())
    return 1;
  if (F.isZero() || G.isZero())
    return 0;
  Variable alpha;
 
#ifdef HAVE_FLINT
  if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G,alpha))
  {
    nmod_poly_t FLINTF, FLINTG;
    convertFacCF2nmod_poly_t (FLINTF, F);
    convertFacCF2nmod_poly_t (FLINTG, G);
    mp_limb_t FLINTresult= nmod_poly_resultant (FLINTF, FLINTG);
    nmod_poly_clear (FLINTF);
    nmod_poly_clear (FLINTG);
    return CanonicalForm ((long) FLINTresult);
  }
  return resultant (F, G, F.mvar());
#elif defined(HAVE_NTL)
  if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G,alpha))
  {
    if (fac_NTL_char != getCharacteristic())
    {
      fac_NTL_char= getCharacteristic();
      zz_p::init (getCharacteristic());
    }
    zz_pX NTLF= convertFacCF2NTLzzpX (F);
    zz_pX NTLG= convertFacCF2NTLzzpX (G);
 
    zz_p NTLResult= resultant (NTLF, NTLG);
 
    return CanonicalForm (to_long (rep (NTLResult)));
  }
  //at this point F or G has an algebraic var.
  if (fac_NTL_char != getCharacteristic())
  {
    fac_NTL_char= getCharacteristic();
    zz_p::init (getCharacteristic());
  }
  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
  zz_pE::init (NTLMipo);
  zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
  zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
  zz_pE NTLResult= resultant (NTLF, NTLG);
 
  return convertNTLzzpE2CF (NTLResult, alpha);
#else
  return resultant (F, G, F.mvar());
#endif
}

Variable Documentation

both_non_zero

int both_non_zero = 0

Definition at line 58 of file cfModResultant.cc.

both_zero

int both_zero = 0

Definition at line 61 of file cfModResultant.cc.

degsf

degsf = NEW_ARRAY(int,n + 1)

Definition at line 49 of file cfModResultant.cc.

degsfx

int degsfx

Definition at line 62 of file cfModResultant.cc.

degsg

degsg = NEW_ARRAY(int,n + 1)

Definition at line 50 of file cfModResultant.cc.

degsgx

int degsgx

Definition at line 62 of file cfModResultant.cc.

else

else

Initial value:

{
    
    for (int i= 1; i <= n; i++)
    {
      if (degsf[i] == 0 && degsg[i] == 0)
      {
        both_zero++;
        continue;
      }
      else
      {
        if (both_zero != 0 && i != 1)
        {
          M.newpair (Variable (i), Variable (i - both_zero));
          N.newpair (Variable (i - both_zero), Variable (i));
        }
      }
    }
  }
 
  DELETE_ARRAY(degsf)

Definition at line 216 of file cfModResultant.cc.

f_zero

int f_zero = 0

Definition at line 59 of file cfModResultant.cc.

Flevel

int Flevel =F.level()

Definition at line 63 of file cfModResultant.cc.

G

const CanonicalForm& G

Definition at line 45 of file cfModResultant.cc.

g_zero

int g_zero = 0

Definition at line 60 of file cfModResultant.cc.

Glevel

int Glevel =G.level()

Definition at line 64 of file cfModResultant.cc.

M

const CanonicalForm CFMap& M

Definition at line 45 of file cfModResultant.cc.

N

const CanonicalForm CFMap CFMap& N

Definition at line 46 of file cfModResultant.cc.

x

const CanonicalForm CFMap CFMap const Variable& x

Initial value:

{
  int n= tmax (F.level(), G.level())

Definition at line 46 of file cfModResultant.cc.