My Project
Loading...
Searching...
No Matches
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
22gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
23{
25 ZZX G1=convertFacCF2NTLZZX(G);
26 ZZX R=GCD(F1,G1);
27 return convertNTLZZX2CF(R,F.mvar());
28}
29
31gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
32{
34 {
36 zz_p::init(getCharacteristic());
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);
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
75CanonicalForm 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 }
116 Dp = gcd( mapinto( f ), mapinto( g ) );
117 Dp = ( Dp / Dp.lc() ) * mapinto( cg );
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());
203 nmod_poly_clear (G1);
207 return r;
208 }
209#elif defined(HAVE_NTL)
211 && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
212 {
214 {
216 zz_p::init(getCharacteristic());
217 }
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);
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 FACTORY_PUBLIC content(const CanonicalForm &)
CanonicalForm content ( const CanonicalForm & f )
Definition: cf_gcd.cc:603
CanonicalForm mapinto(const CanonicalForm &f)
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
CanonicalForm FACTORY_PUBLIC pp(const CanonicalForm &)
CanonicalForm pp ( const CanonicalForm & f )
Definition: cf_gcd.cc:676
void FACTORY_PUBLIC setCharacteristic(int c)
Definition: cf_char.cc:28
CanonicalForm FACTORY_PUBLIC gcd_poly(const CanonicalForm &f, const CanonicalForm &g)
CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
Definition: cf_gcd.cc:492
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:4078
cl
Definition: cfModGcd.cc:4100
g
Definition: cfModGcd.cc:4090
CanonicalForm b
Definition: cfModGcd.cc:4103
CanonicalForm cg
Definition: cfModGcd.cc:4083
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:18
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
#define R
Definition: sirandom.c:27
#define A
Definition: sirandom.c:24
#define M
Definition: sirandom.c:25
int F1(int a1, int &r1)
F1.
int gcd(int a, int b)
Definition: walkSupport.cc:836