My Project  debian-1:4.1.1-p2+ds-4build3
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 
10 #ifdef HAVE_NTL
11 #include "NTLconvert.h"
12 #endif
13 
14 #ifdef HAVE_FLINT
15 #include "FLINTconvert.h"
16 #endif
17 
18 #ifdef HAVE_NTL
19 #ifndef HAVE_FLINT
21 gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
22 {
23  ZZX F1=convertFacCF2NTLZZX(F);
24  ZZX G1=convertFacCF2NTLZZX(G);
25  ZZX R=GCD(F1,G1);
26  return convertNTLZZX2CF(R,F.mvar());
27 }
28 
30 gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
31 {
33  {
35  zz_p::init(getCharacteristic());
36  }
37  zz_pX F1=convertFacCF2NTLzzpX(F);
38  zz_pX G1=convertFacCF2NTLzzpX(G);
39  zz_pX R=GCD(F1,G1);
40  return convertNTLzzpX2CF(R,F.mvar());
41 }
42 #endif
43 #endif
44 
45 #ifdef HAVE_FLINT
48 {
49  nmod_poly_t F1, G1;
52  nmod_poly_gcd (F1, F1, G1);
54  nmod_poly_clear (F1);
55  nmod_poly_clear (G1);
56  return result;
57 }
58 
61 {
62  fmpz_poly_t F1, G1;
65  fmpz_poly_gcd (F1, F1, G1);
67  fmpz_poly_clear (F1);
68  fmpz_poly_clear (G1);
69  return result;
70 }
71 #endif
72 
73 #ifndef HAVE_NTL
74 CanonicalForm gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
75 {
76  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
77  int p, i;
78 
79  if ( primitive )
80  {
81  f = F;
82  g = G;
83  c = 1;
84  }
85  else
86  {
87  CanonicalForm cF = content( F ), cG = content( G );
88  f = F / cF;
89  g = G / cG;
90  c = bgcd( cF, cG );
91  }
92  cg = gcd( f.lc(), g.lc() );
93  cl = ( f.lc() / cg ) * g.lc();
94 // B = 2 * cg * tmin(
95 // maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
96 // maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
97 // )+1;
98  M = tmin( maxNorm(f), maxNorm(g) );
99  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
100  q = 0;
101  i = cf_getNumSmallPrimes() - 1;
102  while ( true )
103  {
104  B = BB;
105  while ( i >= 0 && q < B )
106  {
107  p = cf_getSmallPrime( i );
108  i--;
109  while ( i >= 0 && mod( cl, p ) == 0 )
110  {
111  p = cf_getSmallPrime( i );
112  i--;
113  }
114  setCharacteristic( p );
115  Dp = gcd( mapinto( f ), mapinto( g ) );
116  Dp = ( Dp / Dp.lc() ) * mapinto( cg );
117  setCharacteristic( 0 );
118  if ( Dp.degree() == 0 )
119  return c;
120  if ( q.isZero() )
121  {
122  D = mapinto( Dp );
123  q = p;
124  B = power(CanonicalForm(2),D.degree())*M+1;
125  }
126  else
127  {
128  if ( Dp.degree() == D.degree() )
129  {
130  chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
131  q = newq;
132  D = newD;
133  }
134  else if ( Dp.degree() < D.degree() )
135  {
136  // all previous p's are bad primes
137  q = p;
138  D = mapinto( Dp );
139  B = power(CanonicalForm(2),D.degree())*M+1;
140  }
141  // else p is a bad prime
142  }
143  }
144  if ( i >= 0 )
145  {
146  // now balance D mod q
147  D = pp( balance_p( D, q ) );
148  if ( fdivides( D, f ) && fdivides( D, g ) )
149  return D * c;
150  else
151  q = 0;
152  }
153  else
154  return gcd_poly( F, G );
155  DEBOUTLN( cerr, "another try ..." );
156  }
157 }
158 #endif
159 
160 /** CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
161  *
162  * extgcd() - returns polynomial extended gcd of f and g.
163  *
164  * Returns gcd(f, g) and a and b sucht that f*a+g*b=gcd(f, g).
165  * The gcd is calculated using an extended euclidean polynomial
166  * remainder sequence, so f and g should be polynomials over an
167  * euclidean domain. Normalizes result.
168  *
169  * Note: be sure that f and g have the same level!
170  *
171 **/
174 {
175  if (f.isZero())
176  {
177  a= 0;
178  b= 1;
179  return g;
180  }
181  else if (g.isZero())
182  {
183  a= 1;
184  b= 0;
185  return f;
186  }
187 #ifdef HAVE_NTL
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 #else
211  && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
212  {
214  {
216  zz_p::init(getCharacteristic());
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 #else
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 #endif
305  // may contain bug in the co-factors, see track 107
306  CanonicalForm contf = content( f );
307  CanonicalForm contg = content( g );
308 
309  CanonicalForm p0 = f / contf, p1 = g / contg;
310  CanonicalForm f0 = 1, f1 = 0, g0 = 0, g1 = 1, q, r;
311 
312  while ( ! p1.isZero() )
313  {
314  divrem( p0, p1, q, r );
315  p0 = p1; p1 = r;
316  r = g0 - g1 * q;
317  g0 = g1; g1 = r;
318  r = f0 - f1 * q;
319  f0 = f1; f1 = r;
320  }
321  CanonicalForm contp0 = content( p0 );
322  a = f0 / ( contf * contp0 );
323  b = g0 / ( contg * contp0 );
324  p0 /= contp0;
325  if ( p0.sign() < 0 )
326  {
327  p0 = -p0;
328  a = -a;
329  b = -b;
330  }
331  return p0;
332 }
fac_NTL_char
long fac_NTL_char
Definition: NTLconvert.cc:41
convertZZ2CF
CanonicalForm convertZZ2CF(const ZZ &a)
NAME: convertZZ2CF.
Definition: NTLconvert.cc:489
convertnmod_poly_t2FacCF
CanonicalForm convertnmod_poly_t2FacCF(const nmod_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z/p to CanonicalForm
Definition: FLINTconvert.cc:137
convertFacCF2NTLzzpX
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
Definition: NTLconvert.cc:100
f
FILE * f
Definition: checklibs.c:9
DEBOUTLN
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
result
return result
Definition: facAbsBiFact.cc:76
nmod_poly_clear
nmod_poly_clear(FLINTmipo)
CFFactory::gettype
static int gettype()
Definition: cf_factory.h:28
power
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Definition: canonicalform.cc:1837
g
g
Definition: cfModGcd.cc:4031
cfGcdUtil.h
mod
CF_NO_INLINE CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
Definition: cf_inline.cc:564
CanonicalForm::sign
int sign() const
int CanonicalForm::sign () const
Definition: canonicalform.cc:1295
isPurePoly
bool isPurePoly(const CanonicalForm &f)
Definition: cf_factor.cc:229
cf_primes.h
ftmpl_functions.h
getCharacteristic
int getCharacteristic()
Definition: cf_char.cc:51
b
CanonicalForm b
Definition: cfModGcd.cc:4044
extgcd
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:173
cl
cl
Definition: cfModGcd.cc:4041
CanonicalForm
factory's main class
Definition: canonicalform.h:83
maxNorm
CanonicalForm maxNorm(const CanonicalForm &f)
CanonicalForm maxNorm ( const CanonicalForm & f )
Definition: cf_algorithm.cc:534
divrem
void divrem(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r)
Definition: canonicalform.cc:967
GaloisFieldDomain
#define GaloisFieldDomain
Definition: cf_defs.h:22
i
int i
Definition: cfEzgcd.cc:125
convertFacCF2nmod_poly_t
convertFacCF2nmod_poly_t(FLINTmipo, M)
M
#define M
Definition: sirandom.c:24
content
CanonicalForm content(const CanonicalForm &)
CanonicalForm content ( const CanonicalForm & f )
Definition: cf_gcd.cc:180
nmod_poly_init
nmod_poly_init(FLINTmipo, getCharacteristic())
D
#define D(A)
Definition: gentable.cc:129
gcd_univar_flint0
CanonicalForm gcd_univar_flint0(const CanonicalForm &F, const CanonicalForm &G)
Definition: cfUnivarGcd.cc:60
setCharacteristic
void setCharacteristic(int c)
Definition: cf_char.cc:23
NTLconvert.h
cg
CanonicalForm cg
Definition: cfModGcd.cc:4024
cf_algorithm.h
tmin
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)
chineseRemainder
void 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:52
pp
CanonicalForm pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:253
FLINTconvert.h
convertFacCF2NTLZZX
ZZX convertFacCF2NTLZZX(const CanonicalForm &f)
Definition: NTLconvert.cc:685
fdivides
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
Definition: cf_algorithm.cc:338
convertFacCF2Fmpq_poly_t
void convertFacCF2Fmpq_poly_t(fmpq_poly_t result, const CanonicalForm &f)
conversion of a factory univariate polynomials over Q to fmpq_poly_t
Definition: FLINTconvert.cc:238
convertFmpz_poly_t2FacCF
CanonicalForm convertFmpz_poly_t2FacCF(const fmpz_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z to CanonicalForm
Definition: FLINTconvert.cc:101
bCommonDen
CanonicalForm bCommonDen(const CanonicalForm &f)
CanonicalForm bCommonDen ( const CanonicalForm & f )
Definition: cf_algorithm.cc:293
CanonicalForm::inCoeffDomain
bool inCoeffDomain() const
Definition: canonicalform.cc:119
CanonicalForm::degree
int degree() const
Returns -1 for the zero polynomial and 0 if CO is in a base domain.
Definition: canonicalform.cc:381
convertFmpq_poly_t2FacCF
CanonicalForm convertFmpq_poly_t2FacCF(const fmpq_poly_t p, const Variable &x)
conversion of a FLINT poly over Q to CanonicalForm
Definition: FLINTconvert.cc:212
cf_getSmallPrime
int cf_getSmallPrime(int i)
Definition: cf_primes.cc:28
B
b *CanonicalForm B
Definition: facBivar.cc:52
gcd_univar_flintp
CanonicalForm gcd_univar_flintp(const CanonicalForm &F, const CanonicalForm &G)
Definition: cfUnivarGcd.cc:47
gcd_poly
CanonicalForm gcd_poly(const CanonicalForm &f, const CanonicalForm &g)
CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
Definition: cf_gcd.cc:94
CanonicalForm::mvar
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
Definition: canonicalform.cc:560
mapinto
CanonicalForm mapinto(const CanonicalForm &f)
Definition: canonicalform.h:348
CanonicalForm::lc
CanonicalForm lc() const
CanonicalForm CanonicalForm::lc (), Lc (), LC (), LC ( v ) const.
Definition: canonicalform.cc:304
balance_p
CanonicalForm balance_p(const CanonicalForm &f, const CanonicalForm &q, const CanonicalForm &qh)
same as balance_p ( const CanonicalForm & f, const CanonicalForm & q ) but qh= q/2 is provided,...
Definition: cfGcdUtil.cc:227
bgcd
CanonicalForm bgcd(const CanonicalForm &f, const CanonicalForm &g)
CanonicalForm bgcd ( const CanonicalForm & f, const CanonicalForm & g )
Definition: canonicalform.cc:1589
R
#define R
Definition: sirandom.c:26
gcd
int gcd(int a, int b)
Definition: walkSupport.cc:836
convertNTLZZX2CF
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
Definition: NTLconvert.cc:278
p
int p
Definition: cfModGcd.cc:4019
CanonicalForm::isZero
CF_NO_INLINE bool isZero() const
Definition: cf_inline.cc:372
cf_getNumSmallPrimes
int cf_getNumSmallPrimes()
Definition: cf_primes.cc:34
convertFacCF2Fmpz_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
Definition: FLINTconvert.cc:74
G
static TreeM * G
Definition: janet.cc:32
A
#define A
Definition: sirandom.c:23
convertNTLzzpX2CF
CanonicalForm convertNTLzzpX2CF(const zz_pX &poly, const Variable &x)
Definition: NTLconvert.cc:248
debug.h