My Project
Functions
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]Asome poly
[in]Bsome poly
[in]xsome polynomial variable
[in]probif true use probabilistic algorithm

Definition at line 348 of file cfModResultant.cc.

350 {
351  ASSERT (getCharacteristic() > 0, "characteristic > 0 expected");
352 
353  if (A.isZero() || B.isZero())
354  return 0;
355 
356  int degAx= degree (A, x);
357  int degBx= degree (B, x);
358  if (A.level() < x.level())
359  return power (A, degBx);
360  if (B.level() < x.level())
361  return power (B, degAx);
362 
363  if (degAx == 0)
364  return power (A, degBx);
365  else if (degBx == 0)
366  return power (B, degAx);
367 
368  if (A.isUnivariate() && B.isUnivariate() && A.level() == B.level())
369  return uniResultant (A, B);
370 
371  CanonicalForm F= A;
372  CanonicalForm G= B;
373 
374  CFMap M, N;
375  myCompress (F, G, M, N, x);
376 
377  F= M (F);
378  G= M (G);
379 
380  Variable y= Variable (2);
381 
382  CanonicalForm GEval, FEval, recResult, H;
383  CanonicalForm newtonPoly= 1;
384  CanonicalForm modResult= 0;
385 
386  Variable z= Variable (1);
387  int bound= degAx*degree (G, 2) + degree (F, 2)*degBx;
388 
389  int p= getCharacteristic();
390  CanonicalForm minpoly;
391  Variable alpha= Variable (tmax (F.level(), G.level()) + 1);
392  bool algExt= hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha);
393  CFGenerator * gen;
394  bool extOfExt= false;
395  Variable v= alpha;
396  CanonicalForm primElemAlpha, imPrimElemAlpha;
397  CFList source,dest;
398  if (!algExt && (p < (1 << 28)))
399  {
400  // pass to an extension of size at least 2^29
401  // for very very large input that is maybe too small though
402  int deg= ceil (29.0*((double) log (2)/log (p)))+1;
403  minpoly= randomIrredpoly (deg, z);
404  alpha= rootOf (minpoly);
405  AlgExtGenerator AlgExtGen (alpha);
406  gen= AlgExtGen.clone();
407  for (int i= 0; i < p; i++) // skip values from the prime field
408  (*gen).next();
409  }
410  else if (!algExt)
411  {
412  FFGenerator FFGen;
413  gen= FFGen.clone();
414  }
415  else
416  {
417  int deg= ceil (29.0*((double) log (2)/log (p)));
418  if (degree (getMipo (alpha)) < deg)
419  {
420  mpz_t field_size;
421  mpz_init (field_size);
422  mpz_ui_pow_ui (field_size, p,
423  deg + degree (getMipo (alpha)) - deg%degree (getMipo (alpha)));
424 
425  // field_size needs to fit in an int because of mapUp, mapDown, length of lists etc.
426  if (mpz_fits_sint_p (field_size))
427  {
428  minpoly= randomIrredpoly (deg + degree (getMipo (alpha))
429  - deg%degree (getMipo (alpha)), z);
430  v= rootOf (minpoly);
431  Variable V_buf2;
432  bool primFail= false;
433  extOfExt= true;
434  primElemAlpha= primitiveElement (alpha, V_buf2, primFail);
435  ASSERT (!primFail, "failure in integer factorizer");
436  if (primFail)
437  ; //ERROR
438  else
439  imPrimElemAlpha= mapPrimElem (primElemAlpha, alpha, v);
440  F= mapUp (F, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
441  G= mapUp (G, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
442  }
443  else
444  {
445  deg= deg - deg % degree (getMipo (alpha));
446  mpz_ui_pow_ui (field_size, p, deg);
447  while (deg / degree (getMipo (alpha)) >= 2 && !mpz_fits_sint_p (field_size))
448  {
449  deg -= degree (getMipo (alpha));
450  mpz_ui_pow_ui (field_size, p, deg);
451  }
452  if (deg != degree (getMipo (alpha)))
453  {
454  minpoly= randomIrredpoly (deg, z);
455  v= rootOf (minpoly);
456  Variable V_buf2;
457  bool primFail= false;
458  extOfExt= true;
459  primElemAlpha= primitiveElement (alpha, V_buf2, primFail);
460  ASSERT (!primFail, "failure in integer factorizer");
461  if (primFail)
462  ; //ERROR
463  else
464  imPrimElemAlpha= mapPrimElem (primElemAlpha, alpha, v);
465  F= mapUp (F, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
466  G= mapUp (G, alpha, v, primElemAlpha, imPrimElemAlpha, source, dest);
467  }
468  }
469  mpz_clear (field_size);
470  }
471  AlgExtGenerator AlgExtGen (v);
472  gen= AlgExtGen.clone();
473  for (int i= 0; i < p; i++)
474  (*gen).next();
475  }
476  int count= 0;
477  int equalCount= 0;
478  CanonicalForm point;
479  do
480  {
481  evalPoint (F, G, FEval, GEval, *gen);
482 
483  recResult= resultantFp (FEval, GEval, z, prob);
484 
485  H= newtonInterp ((*gen).item(), recResult, newtonPoly, modResult, y);
486 
487  if (H == modResult)
488  equalCount++;
489  else
490  equalCount= 0;
491 
492  count++;
493  if (count > bound || (prob && equalCount == 2 && !H.inCoeffDomain()))
494  {
495  if (!algExt && degree (H, alpha) <= 0)
496  break;
497  else if (algExt)
498  {
499  if (extOfExt && !isInExtension (H, imPrimElemAlpha, 1, primElemAlpha,
500  dest, source))
501  {
502  H= mapDown (H, primElemAlpha, imPrimElemAlpha, alpha, dest, source);
503  prune (v);
504  break;
505  }
506  else if (!extOfExt)
507  break;
508  }
509  }
510 
511  modResult= H;
512  newtonPoly *= (y - (*gen).item());
513  if ((*gen).hasItems())
514  (*gen).next();
515  else
516  STICKYASSERT (0, "out of evaluation points");
517  } while (1);
518 
519  delete gen;
520 
521  return N (H);
522 }
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
int degree(const CanonicalForm &f)
bool hasFirstAlgVar(const CanonicalForm &f, Variable &a)
check if poly f contains an algebraic variable a
Definition: cf_ops.cc:679
int FACTORY_PUBLIC getCharacteristic()
Definition: cf_char.cc:70
int i
Definition: cfEzgcd.cc:132
int myCompress(const CanonicalForm &F, const CanonicalForm &G, CFMap &M, CFMap &N, bool topLevel)
compressing two polynomials F and G, M is used for compressing, N to reverse the compression
Definition: cfModGcd.cc:93
int p
Definition: cfModGcd.cc:4080
const CanonicalForm CFMap CFMap & N
const CanonicalForm CFMap CFMap const Variable & x
CanonicalForm resultantFp(const CanonicalForm &A, const CanonicalForm &B, const Variable &x, bool prob)
modular resultant algorihtm over Fp
static CanonicalForm uniResultant(const CanonicalForm &F, const CanonicalForm &G)
static void evalPoint(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &FEval, CanonicalForm &GEval, CFGenerator &evalPoint)
const CanonicalForm & G
const CanonicalForm CFMap & M
static CanonicalForm newtonInterp(const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x)
#define STICKYASSERT(expression, message)
Definition: cf_assert.h:64
#define ASSERT(expression, message)
Definition: cf_assert.h:99
CanonicalForm randomIrredpoly(int i, const Variable &x)
computes a random monic irreducible univariate polynomial in x over Fp of degree i via NTL/FLINT
Definition: cf_irred.cc:26
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
CanonicalForm mapPrimElem(const CanonicalForm &primElem, const Variable &alpha, const Variable &beta)
compute the image of a primitive element of in . We assume .
Definition: cf_map_ext.cc:450
CanonicalForm primitiveElement(const Variable &alpha, Variable &beta, bool &fail)
determine a primitive element of , is a primitive element of a field which is isomorphic to
Definition: cf_map_ext.cc:342
static CanonicalForm mapDown(const CanonicalForm &F, const Variable &alpha, const CanonicalForm &G, CFList &source, CFList &dest)
the CanonicalForm G is the output of map_up, returns F considered as an element over ,...
Definition: cf_map_ext.cc:123
static CanonicalForm mapUp(const Variable &alpha, const Variable &beta)
and is a primitive element, returns the image of
Definition: cf_map_ext.cc:70
generate all elements in F_p(alpha) starting from 0
Definition: cf_generator.h:94
virtual class for generators
Definition: cf_generator.h:22
virtual CFGenerator * clone() const
Definition: cf_generator.h:30
virtual void next()
Definition: cf_generator.h:29
class CFMap
Definition: cf_map.h:85
factory's main class
Definition: canonicalform.h:86
CF_NO_INLINE bool isZero() const
bool inCoeffDomain() const
int level() const
level() returns the level of CO.
bool isUnivariate() const
generate all elements in F_p starting from 0
Definition: cf_generator.h:56
CFGenerator * clone() const
Definition: cf_generator.cc:52
factory's class for variables
Definition: factory.h:134
Variable alpha
Definition: facAbsBiFact.cc:51
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53
CanonicalForm H
Definition: facAbsFact.cc:60
b *CanonicalForm B
Definition: facBivar.cc:52
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
bool isInExtension(const CanonicalForm &F, const CanonicalForm &gamma, const int k, const CanonicalForm &delta, CFList &source, CFList &dest)
tests if F is not contained in a subfield defined by gamma (Fq case) or k (GF case)
Variable FACTORY_PUBLIC rootOf(const CanonicalForm &, char name='@')
returns a symbolic root of polynomial with name name Use it to define algebraic variables
Definition: variable.cc:162
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition: variable.cc:207
void FACTORY_PUBLIC prune(Variable &alpha)
Definition: variable.cc:261
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
gmp_float log(const gmp_float &a)
Definition: mpr_complex.cc:343
const signed long ceil(const ampf< Precision > &x)
Definition: amp.h:788
int status int void size_t count
Definition: si_signals.h:59
#define A
Definition: sirandom.c:24

◆ 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]Asome poly
[in]Bsome poly
[in]xsome polynomial variable
[in]probif true use probabilistic algorithm

Definition at line 559 of file cfModResultant.cc.

561 {
562  ASSERT (getCharacteristic() == 0, "characteristic > 0 expected");
563 #ifndef NOASSERT
564  bool isRat= isOn (SW_RATIONAL);
565  On (SW_RATIONAL);
566  ASSERT (bCommonDen (A).isOne(), "input A is rational");
567  ASSERT (bCommonDen (B).isOne(), "input B is rational");
568  if (!isRat)
569  Off (SW_RATIONAL);
570 #endif
571 
572  int degAx= degree (A, x);
573  int degBx= degree (B, x);
574  if (A.level() < x.level())
575  return power (A, degBx);
576  if (B.level() < x.level())
577  return power (B, degAx);
578 
579  if (degAx == 0)
580  return power (A, degBx);
581  else if (degBx == 0)
582  return power (B, degAx);
583 
584  CanonicalForm F= A;
585  CanonicalForm G= B;
586 
587  Variable X= x;
588  if (F.level() != x.level() || G.level() != x.level())
589  {
590  if (F.level() > G.level())
591  X= F.mvar();
592  else
593  X= G.mvar();
594  F= swapvar (F, X, x);
595  G= swapvar (G, X, x);
596  }
597 
598  // now X is the main variable
599 
600  CanonicalForm d= 0;
601  CanonicalForm dd= 0;
603  for (CFIterator i= F; i.hasTerms(); i++)
604  {
605  buf= oneNorm (i.coeff());
606  d= (buf > d) ? buf : d;
607  }
608  CanonicalForm e= 0, ee= 0;
609  for (CFIterator i= G; i.hasTerms(); i++)
610  {
611  buf= oneNorm (i.coeff());
612  e= (buf > e) ? buf : e;
613  }
614  d= power (d, degBx);
615  e= power (e, degAx);
616  CanonicalForm bound= 1;
617  for (int i= degBx + degAx; i > 1; i--)
618  bound *= i;
619  bound *= d*e;
620  bound *= 2;
621 
622  bool onRational= isOn (SW_RATIONAL);
623  if (onRational)
624  Off (SW_RATIONAL);
625  int i = cf_getNumBigPrimes() - 1;
626  int p;
627  CanonicalForm l= lc (F)*lc(G);
628  CanonicalForm resultModP, q (0), newResult, newQ;
630  int equalCount= 0;
631  CanonicalForm test, newTest;
632  int count= 0;
633  do
634  {
635  p = cf_getBigPrime( i );
636  i--;
637  while ( i >= 0 && mod( l, p ) == 0)
638  {
639  p = cf_getBigPrime( i );
640  i--;
641  }
642 
643  if (i <= 0)
644  return resultant (A, B, x);
645 
647 
648  TIMING_START (fac_resultant_p);
649  resultModP= resultantFp (mapinto (F), mapinto (G), X, prob);
650  TIMING_END_AND_PRINT (fac_resultant_p, "time to compute resultant mod p: ");
651 
652  setCharacteristic (0);
653 
654  count++;
655  if ( q.isZero() )
656  {
657  result= mapinto(resultModP);
658  q= p;
659  }
660  else
661  {
662  chineseRemainder( result, q, mapinto (resultModP), p, newResult, newQ );
663  q= newQ;
664  result= newResult;
666  if (test != newTest)
667  {
668  newTest= test;
669  equalCount= 0;
670  }
671  else
672  equalCount++;
673  if (newQ > bound || (prob && equalCount == 2))
674  {
675  result= test;
676  break;
677  }
678  }
679  } while (1);
680 
681  if (onRational)
682  On (SW_RATIONAL);
683  return swapvar (result, X, x);
684 }
bool isOn(int sw)
switches
void On(int sw)
switches
void Off(int sw)
switches
CanonicalForm mapinto(const CanonicalForm &f)
CanonicalForm lc(const CanonicalForm &f)
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
void FACTORY_PUBLIC setCharacteristic(int c)
Definition: cf_char.cc:28
CanonicalForm swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition: cf_ops.cc:168
int l
Definition: cfEzgcd.cc:100
CanonicalForm test
Definition: cfModGcd.cc:4098
static CanonicalForm oneNorm(const CanonicalForm &F)
static CanonicalForm symmetricRemainder(const CanonicalForm &f, const CanonicalForm &q)
CanonicalForm bCommonDen(const CanonicalForm &f)
CanonicalForm bCommonDen ( const CanonicalForm & f )
void FACTORY_PUBLIC chineseRemainder(const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2,...
Definition: cf_chinese.cc:57
CanonicalForm FACTORY_PUBLIC resultant(const CanonicalForm &f, const CanonicalForm &g, const Variable &x)
CanonicalForm resultant ( const CanonicalForm & f, const CanonicalForm & g, const Variable & x )
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:30
int cf_getBigPrime(int i)
Definition: cf_primes.cc:39
int cf_getNumBigPrimes()
Definition: cf_primes.cc:45
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
return result
Definition: facAbsBiFact.cc:75
TIMING_END_AND_PRINT(fac_alg_resultant, "time to compute resultant0: ")
TIMING_START(fac_alg_resultant)
int status int void * buf
Definition: si_signals.h:59