My Project
cfUnivarGcd.cc
Go to the documentation of this file.
1 #include "config.h"
2 
3 #include "debug.h"
4 
5 #include "cf_algorithm.h"
7 #include "cf_primes.h"
8 #include "cfGcdUtil.h"
9 #include "cfUnivarGcd.h"
10 
11 #ifdef HAVE_NTL
12 #include "NTLconvert.h"
13 #endif
14 
15 #ifdef HAVE_FLINT
16 #include "FLINTconvert.h"
17 #endif
18 
19 #ifdef HAVE_NTL
20 #ifndef HAVE_FLINT
22 gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
23 {
24  ZZX F1=convertFacCF2NTLZZX(F);
25  ZZX G1=convertFacCF2NTLZZX(G);
26  ZZX R=GCD(F1,G1);
27  return convertNTLZZX2CF(R,F.mvar());
28 }
29 
31 gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
32 {
34  {
37  }
38  zz_pX F1=convertFacCF2NTLzzpX(F);
39  zz_pX G1=convertFacCF2NTLzzpX(G);
40  zz_pX R=GCD(F1,G1);
41  return convertNTLzzpX2CF(R,F.mvar());
42 }
43 #endif
44 #endif
45 
46 #ifdef HAVE_FLINT
49 {
50  nmod_poly_t F1, G1;
53  nmod_poly_gcd (F1, F1, G1);
55  nmod_poly_clear (F1);
56  nmod_poly_clear (G1);
57  return result;
58 }
59 
62 {
63  fmpz_poly_t F1, G1;
66  fmpz_poly_gcd (F1, F1, G1);
68  fmpz_poly_clear (F1);
69  fmpz_poly_clear (G1);
70  return result;
71 }
72 #endif
73 
74 #ifndef HAVE_NTL
75 CanonicalForm gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
76 {
77  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
78  int p, i;
79 
80  if ( primitive )
81  {
82  f = F;
83  g = G;
84  c = 1;
85  }
86  else
87  {
88  CanonicalForm cF = content( F ), cG = content( G );
89  f = F / cF;
90  g = G / cG;
91  c = bgcd( cF, cG );
92  }
93  cg = gcd( f.lc(), g.lc() );
94  cl = ( f.lc() / cg ) * g.lc();
95 // B = 2 * cg * tmin(
96 // maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
97 // maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
98 // )+1;
99  M = tmin( maxNorm(f), maxNorm(g) );
100  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
101  q = 0;
102  i = cf_getNumSmallPrimes() - 1;
103  while ( true )
104  {
105  B = BB;
106  while ( i >= 0 && q < B )
107  {
108  p = cf_getSmallPrime( i );
109  i--;
110  while ( i >= 0 && mod( cl, p ) == 0 )
111  {
112  p = cf_getSmallPrime( i );
113  i--;
114  }
115  setCharacteristic( p );
116  Dp = gcd( mapinto( f ), mapinto( g ) );
117  Dp = ( Dp / Dp.lc() ) * mapinto( cg );
118  setCharacteristic( 0 );
119  if ( Dp.degree() == 0 )
120  return c;
121  if ( q.isZero() )
122  {
123  D = mapinto( Dp );
124  q = p;
125  B = power(CanonicalForm(2),D.degree())*M+1;
126  }
127  else
128  {
129  if ( Dp.degree() == D.degree() )
130  {
131  chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
132  q = newq;
133  D = newD;
134  }
135  else if ( Dp.degree() < D.degree() )
136  {
137  // all previous p's are bad primes
138  q = p;
139  D = mapinto( Dp );
140  B = power(CanonicalForm(2),D.degree())*M+1;
141  }
142  // else p is a bad prime
143  }
144  }
145  if ( i >= 0 )
146  {
147  // now balance D mod q
148  D = pp( balance_p( D, q ) );
149  if ( fdivides( D, f ) && fdivides( D, g ) )
150  return D * c;
151  else
152  q = 0;
153  }
154  else
155  return gcd_poly( F, G );
156  DEBOUTLN( cerr, "another try ..." );
157  }
158 }
159 #endif
160 
161 /** CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
162  *
163  * extgcd() - returns polynomial extended gcd of f and g.
164  *
165  * Returns gcd(f, g) and a and b sucht that f*a+g*b=gcd(f, g).
166  * The gcd is calculated using an extended euclidean polynomial
167  * remainder sequence, so f and g should be polynomials over an
168  * euclidean domain. Normalizes result.
169  *
170  * Note: be sure that f and g have the same level!
171  *
172 **/
175 {
176  if (f.isZero())
177  {
178  a= 0;
179  b= 1;
180  return g;
181  }
182  else if (g.isZero())
183  {
184  a= 1;
185  b= 0;
186  return f;
187  }
188 #ifdef HAVE_FLINT
190  && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
191  {
192  nmod_poly_t F1, G1, A, B, R;
198  nmod_poly_xgcd (R, A, B, F1, G1);
199  a= convertnmod_poly_t2FacCF (A, f.mvar());
200  b= convertnmod_poly_t2FacCF (B, f.mvar());
202  nmod_poly_clear (F1);
203  nmod_poly_clear (G1);
204  nmod_poly_clear (A);
205  nmod_poly_clear (B);
206  nmod_poly_clear (R);
207  return r;
208  }
209 #elif defined(HAVE_NTL)
211  && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
212  {
214  {
217  }
218  zz_pX F1=convertFacCF2NTLzzpX(f);
219  zz_pX G1=convertFacCF2NTLzzpX(g);
220  zz_pX R;
221  zz_pX A,B;
222  XGCD(R,A,B,F1,G1);
223  a=convertNTLzzpX2CF(A,f.mvar());
224  b=convertNTLzzpX2CF(B,f.mvar());
225  return convertNTLzzpX2CF(R,f.mvar());
226  }
227 #endif
228 #ifdef HAVE_FLINT
229  if (( getCharacteristic() ==0) && (f.level()==g.level())
230  && isPurePoly(f) && isPurePoly(g))
231  {
232  fmpq_poly_t F1, G1;
235  fmpq_poly_t R, A, B;
236  fmpq_poly_init (R);
237  fmpq_poly_init (A);
238  fmpq_poly_init (B);
239  fmpq_poly_xgcd (R, A, B, F1, G1);
240  a= convertFmpq_poly_t2FacCF (A, f.mvar());
241  b= convertFmpq_poly_t2FacCF (B, f.mvar());
243  fmpq_poly_clear (F1);
244  fmpq_poly_clear (G1);
245  fmpq_poly_clear (A);
246  fmpq_poly_clear (B);
247  fmpq_poly_clear (R);
248  return r;
249  }
250 #elif defined(HAVE_NTL)
251  if (( getCharacteristic() ==0)
252  && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
253  {
256  ZZX F1=convertFacCF2NTLZZX(f*fc);
257  ZZX G1=convertFacCF2NTLZZX(g*gc);
258  ZZX R=GCD(F1,G1);
259  CanonicalForm r=convertNTLZZX2CF(R,f.mvar());
260  ZZ RR;
261  ZZX A,B;
262  if (r.inCoeffDomain())
263  {
264  XGCD(RR,A,B,F1,G1,1);
266  if(!rr.isZero())
267  {
268  a=convertNTLZZX2CF(A,f.mvar())*fc/rr;
269  b=convertNTLZZX2CF(B,f.mvar())*gc/rr;
270  return CanonicalForm(1);
271  }
272  else
273  {
274  F1 /= R;
275  G1 /= R;
276  XGCD (RR, A,B,F1,G1,1);
277  rr=convertZZ2CF(RR);
278  a=convertNTLZZX2CF(A,f.mvar())*(fc/rr);
279  b=convertNTLZZX2CF(B,f.mvar())*(gc/rr);
280  }
281  }
282  else
283  {
284  XGCD(RR,A,B,F1,G1,1);
286  if (!rr.isZero())
287  {
288  a=convertNTLZZX2CF(A,f.mvar())*fc;
289  b=convertNTLZZX2CF(B,f.mvar())*gc;
290  }
291  else
292  {
293  F1 /= R;
294  G1 /= R;
295  XGCD (RR, A,B,F1,G1,1);
296  rr=convertZZ2CF(RR);
297  a=convertNTLZZX2CF(A,f.mvar())*(fc/rr);
298  b=convertNTLZZX2CF(B,f.mvar())*(gc/rr);
299  }
300  return r;
301  }
302  }
303 #endif
304  // may contain bug in the co-factors, see track 107
305  CanonicalForm contf = content( f );
306  CanonicalForm contg = content( g );
307 
308  CanonicalForm p0 = f / contf, p1 = g / contg;
309  CanonicalForm f0 = 1, f1 = 0, g0 = 0, g1 = 1, q, r;
310 
311  while ( ! p1.isZero() )
312  {
313  divrem( p0, p1, q, r );
314  p0 = p1; p1 = r;
315  r = g0 - g1 * q;
316  g0 = g1; g1 = r;
317  r = f0 - f1 * q;
318  f0 = f1; f1 = r;
319  }
320  CanonicalForm contp0 = content( p0 );
321  a = f0 / ( contf * contp0 );
322  b = g0 / ( contg * contp0 );
323  p0 /= contp0;
324  if ( p0.sign() < 0 )
325  {
326  p0 = -p0;
327  a = -a;
328  b = -b;
329  }
330  return p0;
331 }
CanonicalForm convertFmpq_poly_t2FacCF(const fmpq_poly_t p, const Variable &x)
conversion of a FLINT poly over Q to CanonicalForm
CanonicalForm convertnmod_poly_t2FacCF(const nmod_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z/p to CanonicalForm
CanonicalForm convertFmpz_poly_t2FacCF(const fmpz_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z to CanonicalForm
void convertFacCF2Fmpq_poly_t(fmpq_poly_t result, const CanonicalForm &f)
conversion of a factory univariate polynomials over Q to fmpq_poly_t
void convertFacCF2Fmpz_poly_t(fmpz_poly_t result, const CanonicalForm &f)
conversion of a factory univariate polynomial over Z to a fmpz_poly_t
This file defines functions for conversion to FLINT (www.flintlib.org) and back.
CanonicalForm convertZZ2CF(const ZZ &a)
NAME: convertZZ2CF.
Definition: NTLconvert.cc:495
ZZX convertFacCF2NTLZZX(const CanonicalForm &f)
Definition: NTLconvert.cc:691
CanonicalForm convertNTLzzpX2CF(const zz_pX &poly, const Variable &x)
Definition: NTLconvert.cc:255
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
Definition: NTLconvert.cc:285
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
Definition: NTLconvert.cc:105
VAR long fac_NTL_char
Definition: NTLconvert.cc:46
Conversion to and from NTL.
CanonicalForm bgcd(const CanonicalForm &f, const CanonicalForm &g)
CanonicalForm bgcd ( const CanonicalForm & f, const CanonicalForm & g )
void divrem(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r)
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
CanonicalForm mapinto(const CanonicalForm &f)
CanonicalForm gcd_poly(const CanonicalForm &f, const CanonicalForm &g)
CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
Definition: cf_gcd.cc:492
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
void FACTORY_PUBLIC setCharacteristic(int c)
Definition: cf_char.cc:28
CanonicalForm pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:676
CanonicalForm content(const CanonicalForm &)
CanonicalForm content ( const CanonicalForm & f )
Definition: cf_gcd.cc:603
int FACTORY_PUBLIC getCharacteristic()
Definition: cf_char.cc:70
int i
Definition: cfEzgcd.cc:132
coprimality check and change of representation mod n
int p
Definition: cfModGcd.cc:4080
cl
Definition: cfModGcd.cc:4102
g
Definition: cfModGcd.cc:4092
CanonicalForm b
Definition: cfModGcd.cc:4105
CanonicalForm cg
Definition: cfModGcd.cc:4085
CanonicalForm gcd_univar_flint0(const CanonicalForm &F, const CanonicalForm &G)
Definition: cfUnivarGcd.cc:61
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
CanonicalForm gcd_univar_flintp(const CanonicalForm &F, const CanonicalForm &G)
Definition: cfUnivarGcd.cc:48
univariate Gcd over finite fields and Z, extended GCD over finite fields and Q
CanonicalForm bCommonDen(const CanonicalForm &f)
CanonicalForm bCommonDen ( const CanonicalForm & f )
CanonicalForm maxNorm(const CanonicalForm &f)
CanonicalForm maxNorm ( const CanonicalForm & f )
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
declarations of higher level algorithms.
bool isPurePoly(const CanonicalForm &f)
Definition: cf_factor.cc:244
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
#define GaloisFieldDomain
Definition: cf_defs.h:24
static CanonicalForm gcd_poly_univar0(const CanonicalForm &F, const CanonicalForm &G, bool primitive)
Definition: cf_gcd.cc:178
static CanonicalForm balance_p(const CanonicalForm &f, const CanonicalForm &q, const CanonicalForm &qh)
Definition: cf_gcd.cc:149
int cf_getNumSmallPrimes()
Definition: cf_primes.cc:34
int cf_getSmallPrime(int i)
Definition: cf_primes.cc:28
access to prime tables
FILE * f
Definition: checklibs.c:9
static int gettype()
Definition: cf_factory.h:28
factory's main class
Definition: canonicalform.h:86
CF_NO_INLINE bool isZero() const
int degree() const
Returns -1 for the zero polynomial and 0 if CO is in a base domain.
int sign() const
int CanonicalForm::sign () const
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
bool inCoeffDomain() const
CanonicalForm lc() const
CanonicalForm CanonicalForm::lc (), Lc (), LC (), LC ( v ) const.
functions to print debug output
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
return result
Definition: facAbsBiFact.cc:75
b *CanonicalForm B
Definition: facBivar.cc:52
nmod_poly_init(FLINTmipo, getCharacteristic())
convertFacCF2nmod_poly_t(FLINTmipo, M)
nmod_poly_clear(FLINTmipo)
some useful template functions.
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)
#define D(A)
Definition: gentable.cc:131
STATIC_VAR TreeM * G
Definition: janet.cc:31
void init()
Definition: lintree.cc:864
#define R
Definition: sirandom.c:27
#define A
Definition: sirandom.c:24
#define M
Definition: sirandom.c:25
int gcd(int a, int b)
Definition: walkSupport.cc:836