cfModResultant.h File Reference

modular resultant algorithm as described by G. More...

#include "canonicalform.h"

Go to the source code of this file.

Functions
CanonicalForm	resultantFp (const CanonicalForm &A, const CanonicalForm &B, const Variable &x, bool prob=true)
	modular resultant algorihtm over Fp More...
CanonicalForm	resultantZ (const CanonicalForm &A, const CanonicalForm &B, const Variable &x, bool prob=true)
	modular resultant algorihtm over Z More...

Detailed Description

modular resultant algorithm as described by G.

E. Collins in "The Calculation of multivariate polynomial resultants"

Author

Martin Lee

Definition in file cfModResultant.h.

Function Documentation

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);
}