My Project
Functions | Variables
cfGcdAlgExt.cc File Reference
#include "config.h"
#include <stdio.h>
#include <iostream>
#include "cf_assert.h"
#include "timing.h"
#include "templates/ftmpl_functions.h"
#include "cf_defs.h"
#include "canonicalform.h"
#include "cf_iter.h"
#include "cf_primes.h"
#include "cf_algorithm.h"
#include "cfGcdAlgExt.h"
#include "cfUnivarGcd.h"
#include "cf_map.h"
#include "cf_generator.h"
#include "facMul.h"
#include "cfNTLzzpEXGCD.h"
#include "NTLconvert.h"
#include "FLINTconvert.h"

Go to the source code of this file.

Functions

 TIMING_DEFINE_PRINT (alg_content_p) TIMING_DEFINE_PRINT(alg_content) TIMING_DEFINE_PRINT(alg_compress) TIMING_DEFINE_PRINT(alg_termination) TIMING_DEFINE_PRINT(alg_termination_p) TIMING_DEFINE_PRINT(alg_reconstruction) TIMING_DEFINE_PRINT(alg_newton_p) TIMING_DEFINE_PRINT(alg_recursion_p) TIMING_DEFINE_PRINT(alg_gcd_p) TIMING_DEFINE_PRINT(alg_euclid_p) static int myCompress(const CanonicalForm &F
 compressing two polynomials F and G, M is used for compressing, N to reverse the compression More...
 
 for (int i=0;i<=n;i++) degsf[i]
 
 if (topLevel)
 
 DELETE_ARRAY (degsg)
 
void tryInvert (const CanonicalForm &F, const CanonicalForm &M, CanonicalForm &inv, bool &fail)
 
static CanonicalForm trycontent (const CanonicalForm &f, const Variable &x, const CanonicalForm &M, bool &fail)
 
static CanonicalForm tryvcontent (const CanonicalForm &f, const Variable &x, const CanonicalForm &M, bool &fail)
 
static CanonicalForm trycf_content (const CanonicalForm &f, const CanonicalForm &g, const CanonicalForm &M, bool &fail)
 
static CanonicalForm tryNewtonInterp (const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x, const CanonicalForm &M, bool &fail)
 
static void leadDeg (const CanonicalForm &f, int degs[])
 
void tryBrownGCD (const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M, CanonicalForm &result, bool &fail, bool topLevel)
 modular gcd over F_p[x]/(M) for not necessarily irreducible M. If a zero divisor is encountered fail is set to true. More...
 
static CanonicalForm myicontent (const CanonicalForm &f, const CanonicalForm &c)
 
static CanonicalForm myicontent (const CanonicalForm &f)
 
CanonicalForm QGCD (const CanonicalForm &F, const CanonicalForm &G)
 gcd over Q(a) More...
 
bool isLess (int *a, int *b, int lower, int upper)
 
bool isEqual (int *a, int *b, int lower, int upper)
 
CanonicalForm firstLC (const CanonicalForm &f)
 

Variables

const CanonicalFormG
 
const CanonicalForm CFMapM
 
const CanonicalForm CFMap CFMapN
 
const CanonicalForm CFMap CFMap bool topLevel
 
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 Flevel =F.level()
 
int Glevel =G.level()
 
 else
 
 return
 

Function Documentation

◆ DELETE_ARRAY()

DELETE_ARRAY ( degsg  )

◆ firstLC()

CanonicalForm firstLC ( const CanonicalForm f)

Definition at line 955 of file cfGcdAlgExt.cc.

956 { // returns the leading coefficient (LC) of level <= 1
957  CanonicalForm ret = f;
958  while(ret.level() > 1)
959  ret = LC(ret);
960  return ret;
961 }
CanonicalForm LC(const CanonicalForm &f)
FILE * f
Definition: checklibs.c:9
factory's main class
Definition: canonicalform.h:86
int level() const
level() returns the level of CO.

◆ for()

for ( int  i = 0;i<=n;i++)

Definition at line 72 of file cfEzgcd.cc.

73  {
74  if (degsf[i] != 0 && degsg[i] != 0)
75  {
76  both_non_zero++;
77  continue;
78  }
79  if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
80  {
81  f_zero++;
82  continue;
83  }
84  if (degsg[i] == 0 && degsf[i] && i <= F.level())
85  {
86  g_zero++;
87  continue;
88  }
89  }
int * degsf
Definition: cfEzgcd.cc:59
int f_zero
Definition: cfEzgcd.cc:69
const CanonicalForm & G
Definition: cfEzgcd.cc:55
const CanonicalForm CFMap CFMap int & both_non_zero
Definition: cfEzgcd.cc:57
int i
Definition: cfEzgcd.cc:132
int g_zero
Definition: cfEzgcd.cc:70
int * degsg
Definition: cfEzgcd.cc:60

◆ if()

if ( topLevel  )

Definition at line 75 of file cfGcdAlgExt.cc.

76  {
77  for (int i= 1; i <= n; i++)
78  {
79  if (degsf[i] != 0 && degsg[i] != 0)
80  {
81  both_non_zero++;
82  continue;
83  }
84  if (degsf[i] == 0 && degsg[i] != 0 && i <= Glevel)
85  {
86  f_zero++;
87  continue;
88  }
89  if (degsg[i] == 0 && degsf[i] && i <= Flevel)
90  {
91  g_zero++;
92  continue;
93  }
94  }
95 
96  if (both_non_zero == 0)
97  {
100  return 0;
101  }
102 
103  // map Variables which do not occur in both polynomials to higher levels
104  int k= 1;
105  int l= 1;
106  for (int i= 1; i <= n; i++)
107  {
108  if (degsf[i] != 0 && degsg[i] == 0 && i <= Flevel)
109  {
110  if (k + both_non_zero != i)
111  {
112  M.newpair (Variable (i), Variable (k + both_non_zero));
113  N.newpair (Variable (k + both_non_zero), Variable (i));
114  }
115  k++;
116  }
117  if (degsf[i] == 0 && degsg[i] != 0 && i <= Glevel)
118  {
119  if (l + g_zero + both_non_zero != i)
120  {
121  M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
122  N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
123  }
124  l++;
125  }
126  }
127 
128  // sort Variables x_{i} in increasing order of
129  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
130  int m= tmax (Flevel, Glevel);
131  int min_max_deg;
132  k= both_non_zero;
133  l= 0;
134  int i= 1;
135  while (k > 0)
136  {
137  if (degsf [i] != 0 && degsg [i] != 0)
138  min_max_deg= tmax (degsf[i], degsg[i]);
139  else
140  min_max_deg= 0;
141  while (min_max_deg == 0)
142  {
143  i++;
144  min_max_deg= tmax (degsf[i], degsg[i]);
145  if (degsf [i] != 0 && degsg [i] != 0)
146  min_max_deg= tmax (degsf[i], degsg[i]);
147  else
148  min_max_deg= 0;
149  }
150  for (int j= i + 1; j <= m; j++)
151  {
152  if (tmax (degsf[j],degsg[j]) <= min_max_deg && degsf[j] != 0 && degsg [j] != 0)
153  {
154  min_max_deg= tmax (degsf[j], degsg[j]);
155  l= j;
156  }
157  }
158  if (l != 0)
159  {
160  if (l != k)
161  {
162  M.newpair (Variable (l), Variable(k));
163  N.newpair (Variable (k), Variable(l));
164  degsf[l]= 0;
165  degsg[l]= 0;
166  l= 0;
167  }
168  else
169  {
170  degsf[l]= 0;
171  degsg[l]= 0;
172  l= 0;
173  }
174  }
175  else if (l == 0)
176  {
177  if (i != k)
178  {
179  M.newpair (Variable (i), Variable (k));
180  N.newpair (Variable (k), Variable (i));
181  degsf[i]= 0;
182  degsg[i]= 0;
183  }
184  else
185  {
186  degsf[i]= 0;
187  degsg[i]= 0;
188  }
189  i++;
190  }
191  k--;
192  }
193  }
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int k
Definition: cfEzgcd.cc:99
const CanonicalForm CFMap CFMap & N
Definition: cfGcdAlgExt.cc:56
int both_non_zero
Definition: cfGcdAlgExt.cc:68
int * degsf
Definition: cfGcdAlgExt.cc:59
int f_zero
Definition: cfGcdAlgExt.cc:69
int Flevel
Definition: cfGcdAlgExt.cc:72
int Glevel
Definition: cfGcdAlgExt.cc:73
const CanonicalForm CFMap & M
Definition: cfGcdAlgExt.cc:55
int g_zero
Definition: cfGcdAlgExt.cc:70
int * degsg
Definition: cfGcdAlgExt.cc:60
DELETE_ARRAY(degsg)
factory's class for variables
Definition: factory.h:127
int j
Definition: facHensel.cc:110
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)

◆ isEqual()

bool isEqual ( int *  a,
int *  b,
int  lower,
int  upper 
)

Definition at line 946 of file cfGcdAlgExt.cc.

947 { // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
948  for(int i=lower; i<=upper; i++)
949  if(a[i] != b[i])
950  return false;
951  return true;
952 }
CanonicalForm b
Definition: cfModGcd.cc:4103

◆ isLess()

bool isLess ( int *  a,
int *  b,
int  lower,
int  upper 
)

Definition at line 935 of file cfGcdAlgExt.cc.

936 { // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
937  for(int i=upper; i>=lower; i--)
938  if(a[i] == b[i])
939  continue;
940  else
941  return a[i] < b[i];
942  return true;
943 }

◆ leadDeg()

static void leadDeg ( const CanonicalForm f,
int  degs[] 
)
static

Definition at line 372 of file cfGcdAlgExt.cc.

373 { // leading degree vector w.r.t. lex. monomial order x(i+1) > x(i)
374  // 'a' should point to an array of sufficient size level(f)+1
375  if(f.inCoeffDomain())
376  return;
377  CanonicalForm tmp = f;
378  do
379  {
380  degs[tmp.level()] = tmp.degree();
381  tmp = LC(tmp);
382  }
383  while(!tmp.inCoeffDomain());
384 }
int degree() const
Returns -1 for the zero polynomial and 0 if CO is in a base domain.
bool inCoeffDomain() const

◆ myicontent() [1/2]

static CanonicalForm myicontent ( const CanonicalForm f)
static

Definition at line 721 of file cfGcdAlgExt.cc.

722 {
723 #if defined(HAVE_NTL) || defined(HAVE_FLINT)
724  return myicontent( f, 0 );
725 #else
726  return 1;
727 #endif
728 }
static CanonicalForm myicontent(const CanonicalForm &f, const CanonicalForm &c)
Definition: cfGcdAlgExt.cc:671

◆ myicontent() [2/2]

static CanonicalForm myicontent ( const CanonicalForm f,
const CanonicalForm c 
)
static

Definition at line 671 of file cfGcdAlgExt.cc.

672 {
673 #if defined(HAVE_NTL) || defined(HAVE_FLINT)
674  if (f.isOne() || c.isOne())
675  return 1;
676  if ( f.inBaseDomain() && c.inBaseDomain())
677  {
678  if (c.isZero()) return abs(f);
679  return bgcd( f, c );
680  }
681  else if ( (f.inCoeffDomain() && c.inCoeffDomain()) ||
682  (f.inCoeffDomain() && c.inBaseDomain()) ||
683  (f.inBaseDomain() && c.inCoeffDomain()))
684  {
685  if (c.isZero()) return abs (f);
686 #ifdef HAVE_FLINT
687  fmpz_poly_t FLINTf, FLINTc;
688  convertFacCF2Fmpz_poly_t (FLINTf, f);
689  convertFacCF2Fmpz_poly_t (FLINTc, c);
690  fmpz_poly_gcd (FLINTc, FLINTc, FLINTf);
692  if (f.inCoeffDomain())
693  result= convertFmpz_poly_t2FacCF (FLINTc, f.mvar());
694  else
695  result= convertFmpz_poly_t2FacCF (FLINTc, c.mvar());
696  fmpz_poly_clear (FLINTc);
697  fmpz_poly_clear (FLINTf);
698  return result;
699 #else
700  ZZX NTLf= convertFacCF2NTLZZX (f);
701  ZZX NTLc= convertFacCF2NTLZZX (c);
702  NTLc= GCD (NTLc, NTLf);
703  if (f.inCoeffDomain())
704  return convertNTLZZX2CF(NTLc,f.mvar());
705  else
706  return convertNTLZZX2CF(NTLc,c.mvar());
707 #endif
708  }
709  else
710  {
711  CanonicalForm g = c;
712  for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
713  g = myicontent( i.coeff(), g );
714  return g;
715  }
716 #else
717  return 1;
718 #endif
719 }
CanonicalForm convertFmpz_poly_t2FacCF(const fmpz_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z to CanonicalForm
void convertFacCF2Fmpz_poly_t(fmpz_poly_t result, const CanonicalForm &f)
conversion of a factory univariate polynomial over Z to a fmpz_poly_t
Rational abs(const Rational &a)
Definition: GMPrat.cc:436
ZZX convertFacCF2NTLZZX(const CanonicalForm &f)
Definition: NTLconvert.cc:691
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
Definition: NTLconvert.cc:285
CanonicalForm bgcd(const CanonicalForm &f, const CanonicalForm &g)
CanonicalForm bgcd ( const CanonicalForm & f, const CanonicalForm & g )
g
Definition: cfModGcd.cc:4090
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
CF_NO_INLINE bool isZero() const
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
CF_NO_INLINE bool isOne() const
bool inBaseDomain() const
return result
Definition: facAbsBiFact.cc:75

◆ QGCD()

gcd over Q(a)

Definition at line 730 of file cfGcdAlgExt.cc.

731 { // f,g in Q(a)[x1,...,xn]
732  if(F.isZero())
733  {
734  if(G.isZero())
735  return G; // G is zero
736  if(G.inCoeffDomain())
737  return CanonicalForm(1);
738  CanonicalForm lcinv= 1/Lc (G);
739  return G*lcinv; // return monic G
740  }
741  if(G.isZero()) // F is non-zero
742  {
743  if(F.inCoeffDomain())
744  return CanonicalForm(1);
745  CanonicalForm lcinv= 1/Lc (F);
746  return F*lcinv; // return monic F
747  }
748  if(F.inCoeffDomain() || G.inCoeffDomain())
749  return CanonicalForm(1);
750  // here: both NOT inCoeffDomain
751  CanonicalForm f, g, tmp, M, q, D, Dp, cl, newq, mipo;
752  int p, i;
753  int *bound, *other; // degree vectors
754  bool fail;
755  bool off_rational=!isOn(SW_RATIONAL);
756  On( SW_RATIONAL ); // needed by bCommonDen
757  f = F * bCommonDen(F);
758  g = G * bCommonDen(G);
760  CanonicalForm contf= myicontent (f);
761  CanonicalForm contg= myicontent (g);
762  f /= contf;
763  g /= contg;
764  CanonicalForm gcdcfcg= myicontent (contf, contg);
765  TIMING_END_AND_PRINT (alg_content, "time for content in alg gcd: ")
766  Variable a, b;
767  if(hasFirstAlgVar(f,a))
768  {
769  if(hasFirstAlgVar(g,b))
770  {
771  if(b.level() > a.level())
772  a = b;
773  }
774  }
775  else
776  {
777  if(!hasFirstAlgVar(g,a))// both not in extension
778  {
779  Off( SW_RATIONAL );
780  Off( SW_USE_QGCD );
781  tmp = gcdcfcg*gcd( f, g );
782  On( SW_USE_QGCD );
783  if (off_rational) Off(SW_RATIONAL);
784  return tmp;
785  }
786  }
787  // here: a is the biggest alg. var in f and g AND some of f,g is in extension
788  setReduce(a,false); // do not reduce expressions modulo mipo
789  tmp = getMipo(a);
790  M = tmp * bCommonDen(tmp);
791  // here: f, g in Z[a][x1,...,xn], M in Z[a] not necessarily monic
792  Off( SW_RATIONAL ); // needed by mod
793  // calculate upper bound for degree vector of gcd
794  int mv = f.level(); i = g.level();
795  if(i > mv)
796  mv = i;
797  // here: mv is level of the largest variable in f, g
798  bound = new int[mv+1]; // 'bound' could be indexed from 0 to mv, but we will only use from 1 to mv
799  other = new int[mv+1];
800  for(int i=1; i<=mv; i++) // initialize 'bound', 'other' with zeros
801  bound[i] = other[i] = 0;
802  leadDeg(f,bound); // 'bound' is set the leading degree vector of f
803  leadDeg(g,other);
804  for(int i=1; i<=mv; i++) // entry at i=0 not visited
805  if(other[i] < bound[i])
806  bound[i] = other[i];
807  // now 'bound' is the smaller vector
808  cl = lc(M) * lc(f) * lc(g);
809  q = 1;
810  D = 0;
811  CanonicalForm test= 0;
812  bool equal= false;
813  for( i=cf_getNumBigPrimes()-1; i>-1; i-- )
814  {
815  p = cf_getBigPrime(i);
816  if( mod( cl, p ).isZero() ) // skip lc-bad primes
817  continue;
818  fail = false;
820  mipo = mapinto(M);
821  mipo /= mipo.lc();
822  // here: mipo is monic
823  TIMING_START (alg_gcd_p)
824  tryBrownGCD( mapinto(f), mapinto(g), mipo, Dp, fail );
825  TIMING_END_AND_PRINT (alg_gcd_p, "time for alg gcd mod p: ")
826  if( fail ) // mipo splits in char p
827  {
829  continue;
830  }
831  if( Dp.inCoeffDomain() ) // early termination
832  {
833  tryInvert(Dp,mipo,tmp,fail); // check if zero divisor
835  if(fail)
836  continue;
837  setReduce(a,true);
838  if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
839  delete[] other;
840  delete[] bound;
841  return gcdcfcg;
842  }
844  // here: Dp NOT inCoeffDomain
845  for(int i=1; i<=mv; i++)
846  other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
847  leadDeg(Dp,other);
848 
849  if(isEqual(bound, other, 1, mv)) // equal
850  {
851  chineseRemainder( D, q, mapinto(Dp), p, tmp, newq );
852  // tmp = Dp mod p
853  // tmp = D mod q
854  // newq = p*q
855  q = newq;
856  if( D != tmp )
857  D = tmp;
858  On( SW_RATIONAL );
859  TIMING_START (alg_reconstruction)
860  tmp = Farey( D, q ); // Farey
861  tmp *= bCommonDen (tmp);
862  TIMING_END_AND_PRINT (alg_reconstruction,
863  "time for rational reconstruction in alg gcd: ")
864  setReduce(a,true); // reduce expressions modulo mipo
865  On( SW_RATIONAL ); // needed by fdivides
866  if (test != tmp)
867  test= tmp;
868  else
869  equal= true; // modular image did not add any new information
870  TIMING_START (alg_termination)
871 #ifdef HAVE_NTL
872 #ifdef HAVE_FLINT
873  if (equal && tmp.isUnivariate() && f.isUnivariate() && g.isUnivariate()
874  && f.level() == tmp.level() && tmp.level() == g.level())
875  {
876  CanonicalForm Q, R;
877  newtonDivrem (f, tmp, Q, R);
878  if (R.isZero())
879  {
880  newtonDivrem (g, tmp, Q, R);
881  if (R.isZero())
882  {
883  Off (SW_RATIONAL);
884  setReduce (a,true);
885  if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
886  TIMING_END_AND_PRINT (alg_termination,
887  "time for successful termination test in alg gcd: ")
888  delete[] other;
889  delete[] bound;
890  return tmp*gcdcfcg;
891  }
892  }
893  }
894  else
895 #endif
896 #endif
897  if(equal && fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
898  {
899  Off( SW_RATIONAL );
900  setReduce(a,true);
901  if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
902  TIMING_END_AND_PRINT (alg_termination,
903  "time for successful termination test in alg gcd: ")
904  delete[] other;
905  delete[] bound;
906  return tmp*gcdcfcg;
907  }
908  TIMING_END_AND_PRINT (alg_termination,
909  "time for unsuccessful termination test in alg gcd: ")
910  Off( SW_RATIONAL );
911  setReduce(a,false); // do not reduce expressions modulo mipo
912  continue;
913  }
914  if( isLess(bound, other, 1, mv) ) // current prime unlucky
915  continue;
916  // here: isLess(other, bound, 1, mv) ) ==> all previous primes unlucky
917  q = p;
918  D = mapinto(Dp); // shortcut CRA // shortcut CRA
919  for(int i=1; i<=mv; i++) // tighten bound
920  bound[i] = other[i];
921  }
922  // hopefully, we never reach this point
923  setReduce(a,true);
924  delete[] other;
925  delete[] bound;
926  Off( SW_USE_QGCD );
927  D = gcdcfcg*gcd( f, g );
928  On( SW_USE_QGCD );
929  if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
930  return D;
931 }
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
bool hasFirstAlgVar(const CanonicalForm &f, Variable &a)
check if poly f contains an algebraic variable a
Definition: cf_ops.cc:679
CanonicalForm Lc(const CanonicalForm &f)
else
Definition: cfGcdAlgExt.cc:195
bool isLess(int *a, int *b, int lower, int upper)
Definition: cfGcdAlgExt.cc:935
static void leadDeg(const CanonicalForm &f, int degs[])
Definition: cfGcdAlgExt.cc:372
bool isEqual(int *a, int *b, int lower, int upper)
Definition: cfGcdAlgExt.cc:946
for(int i=0;i<=n;i++) degsf[i]
Definition: cfEzgcd.cc:72
if(topLevel)
Definition: cfGcdAlgExt.cc:75
return
Definition: cfGcdAlgExt.cc:218
const CanonicalForm & G
Definition: cfGcdAlgExt.cc:55
void tryInvert(const CanonicalForm &F, const CanonicalForm &M, CanonicalForm &inv, bool &fail)
Definition: cfGcdAlgExt.cc:221
int p
Definition: cfModGcd.cc:4078
return false
Definition: cfModGcd.cc:84
cl
Definition: cfModGcd.cc:4100
CanonicalForm gcdcfcg
Definition: cfModGcd.cc:4101
bool equal
Definition: cfModGcd.cc:4126
CanonicalForm test
Definition: cfModGcd.cc:4096
CanonicalForm bCommonDen(const CanonicalForm &f)
CanonicalForm bCommonDen ( const CanonicalForm & f )
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
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
static const int SW_USE_QGCD
set to 1 to use Encarnacion GCD over Q(a)
Definition: cf_defs.h:43
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:31
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
int cf_getBigPrime(int i)
Definition: cf_primes.cc:39
int cf_getNumBigPrimes()
Definition: cf_primes.cc:45
CanonicalForm lc() const
CanonicalForm CanonicalForm::lc (), Lc (), LC (), LC ( v ) const.
int level() const
Definition: factory.h:143
CanonicalForm mipo
Definition: facAlgExt.cc:57
TIMING_END_AND_PRINT(fac_alg_resultant, "time to compute resultant0: ")
TIMING_START(fac_alg_resultant)
CanonicalForm alg_content(const CanonicalForm &f, const CFList &as)
Definition: facAlgFunc.cc:42
void newtonDivrem(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Q, CanonicalForm &R)
division with remainder of univariate polynomials over Q and Q(a) using Newton inversion,...
Definition: facMul.cc:346
bool isZero(const CFArray &A)
checks if entries of A are zero
void setReduce(const Variable &alpha, bool reduce)
Definition: variable.cc:238
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition: variable.cc:207
static number Farey(number, number, const coeffs)
Definition: flintcf_Q.cc:445
#define D(A)
Definition: gentable.cc:131
STATIC_VAR jList * Q
Definition: janet.cc:30
#define R
Definition: sirandom.c:27
int gcd(int a, int b)
Definition: walkSupport.cc:836

◆ TIMING_DEFINE_PRINT()

TIMING_DEFINE_PRINT ( alg_content_p  ) const &

compressing two polynomials F and G, M is used for compressing, N to reverse the compression

◆ tryBrownGCD()

void tryBrownGCD ( const CanonicalForm F,
const CanonicalForm G,
const CanonicalForm M,
CanonicalForm result,
bool &  fail,
bool  topLevel 
)

modular gcd over F_p[x]/(M) for not necessarily irreducible M. If a zero divisor is encountered fail is set to true.

Definition at line 386 of file cfGcdAlgExt.cc.

387 { // assume F,G are multivariate polys over Z/p(a) for big prime p, M "univariate" polynomial in an algebraic variable
388  // M is assumed to be monic
389  if(F.isZero())
390  {
391  if(G.isZero())
392  {
393  result = G; // G is zero
394  return;
395  }
396  if(G.inCoeffDomain())
397  {
398  tryInvert(G,M,result,fail);
399  if(fail)
400  return;
401  result = 1;
402  return;
403  }
404  // try to make G monic modulo M
405  CanonicalForm inv;
406  tryInvert(Lc(G),M,inv,fail);
407  if(fail)
408  return;
409  result = inv*G;
410  result= reduce (result, M);
411  return;
412  }
413  if(G.isZero()) // F is non-zero
414  {
415  if(F.inCoeffDomain())
416  {
417  tryInvert(F,M,result,fail);
418  if(fail)
419  return;
420  result = 1;
421  return;
422  }
423  // try to make F monic modulo M
424  CanonicalForm inv;
425  tryInvert(Lc(F),M,inv,fail);
426  if(fail)
427  return;
428  result = inv*F;
429  result= reduce (result, M);
430  return;
431  }
432  // here: F,G both nonzero
433  if(F.inCoeffDomain())
434  {
435  tryInvert(F,M,result,fail);
436  if(fail)
437  return;
438  result = 1;
439  return;
440  }
441  if(G.inCoeffDomain())
442  {
443  tryInvert(G,M,result,fail);
444  if(fail)
445  return;
446  result = 1;
447  return;
448  }
449  TIMING_START (alg_compress)
450  CFMap MM,NN;
451  int lev= myCompress (F, G, MM, NN, topLevel);
452  if (lev == 0)
453  {
454  result= 1;
455  return;
456  }
457  CanonicalForm f=MM(F);
458  CanonicalForm g=MM(G);
459  TIMING_END_AND_PRINT (alg_compress, "time to compress in alg gcd: ")
460  // here: f,g are compressed
461  // compute largest variable in f or g (least one is Variable(1))
462  int mv = f.level();
463  if(g.level() > mv)
464  mv = g.level();
465  // here: mv is level of the largest variable in f, g
466  Variable v1= Variable (1);
467 #ifdef HAVE_NTL
468  Variable v= M.mvar();
469  int ch=getCharacteristic();
470  if (fac_NTL_char != ch)
471  {
472  fac_NTL_char= ch;
473  zz_p::init (ch);
474  }
475  zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
476  zz_pE::init (NTLMipo);
477  zz_pEX NTLResult;
478  zz_pEX NTLF;
479  zz_pEX NTLG;
480 #endif
481  if(mv == 1) // f,g univariate
482  {
483  TIMING_START (alg_euclid_p)
484 #ifdef HAVE_NTL
485  NTLF= convertFacCF2NTLzz_pEX (f, NTLMipo);
486  NTLG= convertFacCF2NTLzz_pEX (g, NTLMipo);
487  tryNTLGCD (NTLResult, NTLF, NTLG, fail);
488  if (fail)
489  return;
490  result= convertNTLzz_pEX2CF (NTLResult, f.mvar(), v);
491 #else
492  tryEuclid(f,g,M,result,fail);
493  if(fail)
494  return;
495 #endif
496  TIMING_END_AND_PRINT (alg_euclid_p, "time for euclidean alg mod p: ")
497  result= NN (reduce (result, M)); // do not forget to map back
498  return;
499  }
500  TIMING_START (alg_content_p)
501  // here: mv > 1
502  CanonicalForm cf = tryvcontent(f, Variable(2), M, fail); // cf is univariate poly in var(1)
503  if(fail)
504  return;
505  CanonicalForm cg = tryvcontent(g, Variable(2), M, fail);
506  if(fail)
507  return;
508  CanonicalForm c;
509 #ifdef HAVE_NTL
510  NTLF= convertFacCF2NTLzz_pEX (cf, NTLMipo);
511  NTLG= convertFacCF2NTLzz_pEX (cg, NTLMipo);
512  tryNTLGCD (NTLResult, NTLF, NTLG, fail);
513  if (fail)
514  return;
515  c= convertNTLzz_pEX2CF (NTLResult, v1, v);
516 #else
517  tryEuclid(cf,cg,M,c,fail);
518  if(fail)
519  return;
520 #endif
521  // f /= cf
522  f.tryDiv (cf, M, fail);
523  if(fail)
524  return;
525  // g /= cg
526  g.tryDiv (cg, M, fail);
527  if(fail)
528  return;
529  TIMING_END_AND_PRINT (alg_content_p, "time for content in alg gcd mod p: ")
530  if(f.inCoeffDomain())
531  {
532  tryInvert(f,M,result,fail);
533  if(fail)
534  return;
535  result = NN(c);
536  return;
537  }
538  if(g.inCoeffDomain())
539  {
540  tryInvert(g,M,result,fail);
541  if(fail)
542  return;
543  result = NN(c);
544  return;
545  }
546  int L[mv+1];
547  int N[mv+1];
548  for(int i=2; i<=mv; i++)
549  L[i] = N[i] = 0;
550  leadDeg(f, L);
551  leadDeg(g, N);
552  CanonicalForm gamma;
553  TIMING_START (alg_euclid_p)
554 #ifdef HAVE_NTL
555  NTLF= convertFacCF2NTLzz_pEX (firstLC (f), NTLMipo);
556  NTLG= convertFacCF2NTLzz_pEX (firstLC (g), NTLMipo);
557  tryNTLGCD (NTLResult, NTLF, NTLG, fail);
558  if (fail)
559  return;
560  gamma= convertNTLzz_pEX2CF (NTLResult, v1, v);
561 #else
562  tryEuclid( firstLC(f), firstLC(g), M, gamma, fail );
563  if(fail)
564  return;
565 #endif
566  TIMING_END_AND_PRINT (alg_euclid_p, "time for gcd of lcs in alg mod p: ")
567  for(int i=2; i<=mv; i++) // entries at i=0,1 not visited
568  if(N[i] < L[i])
569  L[i] = N[i];
570  // L is now upper bound for degrees of gcd
571  int dg_im[mv+1]; // for the degree vector of the image we don't need any entry at i=1
572  for(int i=2; i<=mv; i++)
573  dg_im[i] = 0; // initialize
574  CanonicalForm gamma_image, m=1;
575  CanonicalForm gm=0;
576  CanonicalForm g_image, alpha, gnew;
577  FFGenerator gen = FFGenerator();
578  Variable x= Variable (1);
579  bool divides= true;
580  for(FFGenerator gen = FFGenerator(); gen.hasItems(); gen.next())
581  {
582  alpha = gen.item();
583  gamma_image = reduce(gamma(alpha, x),M); // plug in alpha for var(1)
584  if(gamma_image.isZero()) // skip lc-bad points var(1)-alpha
585  continue;
586  TIMING_START (alg_recursion_p)
587  tryBrownGCD( f(alpha, x), g(alpha, x), M, g_image, fail, false ); // recursive call with one var less
588  TIMING_END_AND_PRINT (alg_recursion_p,
589  "time for recursive calls in alg gcd mod p: ")
590  if(fail)
591  return;
592  g_image = reduce(g_image, M);
593  if(g_image.inCoeffDomain()) // early termination
594  {
595  tryInvert(g_image,M,result,fail);
596  if(fail)
597  return;
598  result = NN(c);
599  return;
600  }
601  for(int i=2; i<=mv; i++)
602  dg_im[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
603  leadDeg(g_image, dg_im);
604  if(isEqual(dg_im, L, 2, mv))
605  {
606  CanonicalForm inv;
607  tryInvert (firstLC (g_image), M, inv, fail);
608  if (fail)
609  return;
610  g_image *= inv;
611  g_image *= gamma_image; // multiply by multiple of image lc(gcd)
612  g_image= reduce (g_image, M);
613  TIMING_START (alg_newton_p)
614  gnew= tryNewtonInterp (alpha, g_image, m, gm, x, M, fail);
615  TIMING_END_AND_PRINT (alg_newton_p,
616  "time for Newton interpolation in alg gcd mod p: ")
617  // gnew = gm mod m
618  // gnew = g_image mod var(1)-alpha
619  // mnew = m * (var(1)-alpha)
620  if(fail)
621  return;
622  m *= (x - alpha);
623  if((firstLC(gnew) == gamma) || (gnew == gm)) // gnew did not change
624  {
625  TIMING_START (alg_termination_p)
626  cf = tryvcontent(gnew, Variable(2), M, fail);
627  if(fail)
628  return;
629  divides = true;
630  g_image= gnew;
631  g_image.tryDiv (cf, M, fail);
632  if(fail)
633  return;
634  divides= tryFdivides (g_image,f, M, fail); // trial division (f)
635  if(fail)
636  return;
637  if(divides)
638  {
639  bool divides2= tryFdivides (g_image,g, M, fail); // trial division (g)
640  if(fail)
641  return;
642  if(divides2)
643  {
644  result = NN(reduce (c*g_image, M));
645  TIMING_END_AND_PRINT (alg_termination_p,
646  "time for successful termination test in alg gcd mod p: ")
647  return;
648  }
649  }
650  TIMING_END_AND_PRINT (alg_termination_p,
651  "time for unsuccessful termination test in alg gcd mod p: ")
652  }
653  gm = gnew;
654  continue;
655  }
656 
657  if(isLess(L, dg_im, 2, mv)) // dg_im > L --> current point unlucky
658  continue;
659 
660  // here: isLess(dg_im, L, 2, mv) --> all previous points were unlucky
661  m = CanonicalForm(1); // reset
662  gm = 0; // reset
663  for(int i=2; i<=mv; i++) // tighten bound
664  L[i] = dg_im[i];
665  }
666  // we are out of evaluation points
667  fail = true;
668 }
zz_pEX convertFacCF2NTLzz_pEX(const CanonicalForm &f, const zz_pX &mipo)
Definition: NTLconvert.cc:1064
CanonicalForm convertNTLzz_pEX2CF(const zz_pEX &f, const Variable &x, const Variable &alpha)
Definition: NTLconvert.cc:1092
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
Definition: NTLconvert.cc:105
VAR long fac_NTL_char
Definition: NTLconvert.cc:46
CanonicalForm reduce(const CanonicalForm &f, const CanonicalForm &M)
polynomials in M.mvar() are considered coefficients M univariate monic polynomial the coefficients of...
Definition: cf_ops.cc:660
int level(const CanonicalForm &f)
int FACTORY_PUBLIC getCharacteristic()
Definition: cf_char.cc:70
static CanonicalForm tryNewtonInterp(const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x, const CanonicalForm &M, bool &fail)
Definition: cfGcdAlgExt.cc:357
const CanonicalForm CFMap CFMap bool topLevel
Definition: cfGcdAlgExt.cc:57
static CanonicalForm tryvcontent(const CanonicalForm &f, const Variable &x, const CanonicalForm &M, bool &fail)
CanonicalForm firstLC(const CanonicalForm &f)
Definition: cfGcdAlgExt.cc:955
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:91
Variable x
Definition: cfModGcd.cc:4082
CanonicalForm cf
Definition: cfModGcd.cc:4083
CanonicalForm cg
Definition: cfModGcd.cc:4083
void tryNTLGCD(zz_pEX &x, const zz_pEX &a, const zz_pEX &b, bool &fail)
compute the GCD x of a and b, fail is set to true if a zero divisor is encountered
bool tryFdivides(const CanonicalForm &f, const CanonicalForm &g, const CanonicalForm &M, bool &fail)
same as fdivides but handles zero divisors in Z_p[t]/(f)[x1,...,xn] for reducible f
class CFMap
Definition: cf_map.h:85
CanonicalForm & tryDiv(const CanonicalForm &, const CanonicalForm &, bool &)
same as divremt but handles zero divisors in case we are in Z_p[x]/(f) where f is not irreducible
generate all elements in F_p starting from 0
Definition: cf_generator.h:56
Variable alpha
Definition: facAbsBiFact.cc:51
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
ListNode * next
Definition: janet.h:31
void init()
Definition: lintree.cc:864

◆ trycf_content()

static CanonicalForm trycf_content ( const CanonicalForm f,
const CanonicalForm g,
const CanonicalForm M,
bool &  fail 
)
static

Definition at line 1066 of file cfGcdAlgExt.cc.

1067 { // as cf_content, but takes care of zero divisors
1068  if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
1069  {
1070  CFIterator i = f;
1071  CanonicalForm tmp = g, result;
1072  while ( i.hasTerms() && ! tmp.isOne() && ! fail )
1073  {
1074  tryBrownGCD( i.coeff(), tmp, M, result, fail );
1075  tmp = result;
1076  i++;
1077  }
1078  return result;
1079  }
1080  return abs( f );
1081 }
void tryBrownGCD(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M, CanonicalForm &result, bool &fail, bool topLevel)
modular gcd over F_p[x]/(M) for not necessarily irreducible M. If a zero divisor is encountered fail ...
Definition: cfGcdAlgExt.cc:386
bool getReduce(const Variable &alpha)
Definition: variable.cc:232

◆ trycontent()

static CanonicalForm trycontent ( const CanonicalForm f,
const Variable x,
const CanonicalForm M,
bool &  fail 
)
static

Definition at line 1035 of file cfGcdAlgExt.cc.

1036 { // as 'content', but takes care of zero divisors
1037  ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
1038  Variable y = f.mvar();
1039  if ( y == x )
1040  return trycf_content( f, 0, M, fail );
1041  if ( y < x )
1042  return f;
1043  return swapvar( trycontent( swapvar( f, y, x ), y, M, fail ), y, x );
1044 }
CanonicalForm FACTORY_PUBLIC swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition: cf_ops.cc:168
static CanonicalForm trycontent(const CanonicalForm &f, const Variable &x, const CanonicalForm &M, bool &fail)
static CanonicalForm trycf_content(const CanonicalForm &f, const CanonicalForm &g, const CanonicalForm &M, bool &fail)
#define ASSERT(expression, message)
Definition: cf_assert.h:99
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53

◆ tryInvert()

void tryInvert ( const CanonicalForm F,
const CanonicalForm M,
CanonicalForm inv,
bool &  fail 
)

Definition at line 221 of file cfGcdAlgExt.cc.

222 { // F, M are required to be "univariate" polynomials in an algebraic variable
223  // we try to invert F modulo M
224  if(F.inBaseDomain())
225  {
226  if(F.isZero())
227  {
228  fail = true;
229  return;
230  }
231  inv = 1/F;
232  return;
233  }
235  Variable a = M.mvar();
236  Variable x = Variable(1);
237  if(!extgcd( replacevar( F, a, x ), replacevar( M, a, x ), inv, b ).isOne())
238  fail = true;
239  else
240  inv = replacevar( inv, x, a ); // change back to alg var
241 }
CanonicalForm FACTORY_PUBLIC replacevar(const CanonicalForm &, const Variable &, const Variable &)
CanonicalForm replacevar ( const CanonicalForm & f, const Variable & x1, const Variable & x2 )
Definition: cf_ops.cc:271
CanonicalForm extgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a,...
Definition: cfUnivarGcd.cc:174

◆ tryNewtonInterp()

static CanonicalForm tryNewtonInterp ( const CanonicalForm alpha,
const CanonicalForm u,
const CanonicalForm newtonPoly,
const CanonicalForm oldInterPoly,
const Variable x,
const CanonicalForm M,
bool &  fail 
)
inlinestatic

Definition at line 357 of file cfGcdAlgExt.cc.

360 {
361  CanonicalForm interPoly;
362 
363  CanonicalForm inv;
364  tryInvert (newtonPoly (alpha, x), M, inv, fail);
365  if (fail)
366  return 0;
367 
368  interPoly= oldInterPoly+reduce ((u - oldInterPoly (alpha, x))*inv*newtonPoly, M);
369  return interPoly;
370 }

◆ tryvcontent()

static CanonicalForm tryvcontent ( const CanonicalForm f,
const Variable x,
const CanonicalForm M,
bool &  fail 
)
static

Definition at line 1047 of file cfGcdAlgExt.cc.

1048 { // as vcontent, but takes care of zero divisors
1049  ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
1050  if ( f.mvar() <= x )
1051  return trycontent( f, x, M, fail );
1052  CFIterator i;
1053  CanonicalForm d = 0, e, ret;
1054  for ( i = f; i.hasTerms() && ! d.isOne() && ! fail; i++ )
1055  {
1056  e = tryvcontent( i.coeff(), x, M, fail );
1057  if(fail)
1058  break;
1059  tryBrownGCD( d, e, M, ret, fail );
1060  d = ret;
1061  }
1062  return d;
1063 }

Variable Documentation

◆ both_non_zero

int both_non_zero = 0

Definition at line 68 of file cfGcdAlgExt.cc.

◆ both_zero

int both_zero = 0

Definition at line 71 of file cfGcdAlgExt.cc.

◆ degsf

degsf = NEW_ARRAY(int,n + 1)

Definition at line 59 of file cfGcdAlgExt.cc.

◆ degsg

degsg = NEW_ARRAY(int,n + 1)

Definition at line 60 of file cfGcdAlgExt.cc.

◆ else

else
Initial value:
{
for (int i= 1; i <= n; i++)
{
if (degsf[i] == 0 && degsg[i] == 0)
{
continue;
}
else
{
if (both_zero != 0)
{
M.newpair (Variable (i), Variable (i - both_zero));
N.newpair (Variable (i - both_zero), Variable (i));
}
}
}
}
int both_zero
Definition: cfGcdAlgExt.cc:71

Definition at line 194 of file cfGcdAlgExt.cc.

◆ f_zero

int f_zero = 0

Definition at line 69 of file cfGcdAlgExt.cc.

◆ Flevel

int Flevel =F.level()

Definition at line 72 of file cfGcdAlgExt.cc.

◆ G

Definition at line 55 of file cfGcdAlgExt.cc.

◆ g_zero

int g_zero = 0

Definition at line 70 of file cfGcdAlgExt.cc.

◆ Glevel

int Glevel =G.level()

Definition at line 73 of file cfGcdAlgExt.cc.

◆ M

Definition at line 55 of file cfGcdAlgExt.cc.

◆ N

Definition at line 56 of file cfGcdAlgExt.cc.

◆ return

return

Definition at line 218 of file cfGcdAlgExt.cc.

◆ topLevel

const CanonicalForm CFMap CFMap bool topLevel
Initial value:
{
int n= tmax (F.level(), G.level())

Definition at line 56 of file cfGcdAlgExt.cc.