My Project
Loading...
Searching...
No Matches
cfGcdAlgExt.cc
Go to the documentation of this file.
1
2#include "config.h"
3
4
5#ifndef NOSTREAMIO
6#ifdef HAVE_CSTDIO
7#include <cstdio>
8#else
9#include <stdio.h>
10#endif
11#ifdef HAVE_IOSTREAM_H
12#include <iostream.h>
13#elif defined(HAVE_IOSTREAM)
14#include <iostream>
15#endif
16#endif
17
18#include "cf_assert.h"
19#include "timing.h"
20
22#include "cf_defs.h"
23#include "canonicalform.h"
24#include "cf_iter.h"
25#include "cf_primes.h"
26#include "cf_algorithm.h"
27#include "cfGcdAlgExt.h"
28#include "cfUnivarGcd.h"
29#include "cf_map.h"
30#include "cf_generator.h"
31#include "facMul.h"
32#include "cfNTLzzpEXGCD.h"
33
34#ifdef HAVE_NTL
35#include "NTLconvert.h"
36#endif
37
38#ifdef HAVE_FLINT
39#include "FLINTconvert.h"
40#endif
41
42TIMING_DEFINE_PRINT(alg_content_p)
44TIMING_DEFINE_PRINT(alg_compress)
45TIMING_DEFINE_PRINT(alg_termination)
46TIMING_DEFINE_PRINT(alg_termination_p)
47TIMING_DEFINE_PRINT(alg_reconstruction)
48TIMING_DEFINE_PRINT(alg_newton_p)
49TIMING_DEFINE_PRINT(alg_recursion_p)
50TIMING_DEFINE_PRINT(alg_gcd_p)
51TIMING_DEFINE_PRINT(alg_euclid_p)
52
53/// compressing two polynomials F and G, M is used for compressing,
54/// N to reverse the compression
56 CFMap & N, bool topLevel)
57{
58 int n= tmax (F.level(), G.level());
59 int * degsf= NEW_ARRAY(int,n + 1);
60 int * degsg= NEW_ARRAY(int,n + 1);
61
62 for (int i = 0; i <= n; i++)
63 degsf[i]= degsg[i]= 0;
64
65 degsf= degrees (F, degsf);
66 degsg= degrees (G, degsg);
67
69 int f_zero= 0;
70 int g_zero= 0;
71 int both_zero= 0;
72 int Flevel=F.level();
73 int Glevel=G.level();
74
76 {
77 for (int i= 1; i <= n; i++)
78 {
79 if (degsf[i] != 0 && degsg[i] != 0)
80 {
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;
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 }
194 else
195 {
196 //arrange Variables such that no gaps occur
197 for (int i= 1; i <= n; i++)
198 {
199 if (degsf[i] == 0 && degsg[i] == 0)
200 {
201 both_zero++;
202 continue;
203 }
204 else
205 {
206 if (both_zero != 0)
207 {
208 M.newpair (Variable (i), Variable (i - both_zero));
209 N.newpair (Variable (i - both_zero), Variable (i));
210 }
211 }
212 }
213 }
214
217
218 return 1;
219}
220
221void tryInvert( const CanonicalForm & F, const CanonicalForm & M, CanonicalForm & inv, bool & fail )
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}
242
243#ifndef HAVE_NTL
244void tryDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
246 bool& fail)
247{
248 if (F.inCoeffDomain())
249 {
250 Q= 0;
251 R= F;
252 return;
253 }
254
256 Variable x= F.mvar();
257 A= F;
258 B= G;
259 int degA= degree (A, x);
260 int degB= degree (B, x);
261
262 if (degA < degB)
263 {
264 R= A;
265 Q= 0;
266 return;
267 }
268
269 tryInvert (Lc (B), mipo, inv, fail);
270 if (fail)
271 return;
272
273 R= A;
274 Q= 0;
275 CanonicalForm Qi;
276 for (int i= degA -degB; i >= 0; i--)
277 {
278 if (degree (R, x) == i + degB)
279 {
280 Qi= Lc (R)*inv*power (x, i);
281 Qi= reduce (Qi, mipo);
282 R -= Qi*B;
283 R= reduce (R, mipo);
284 Q += Qi;
285 }
286 }
287}
288
289void tryEuclid( const CanonicalForm & A, const CanonicalForm & B, const CanonicalForm & M, CanonicalForm & result, bool & fail )
290{
292 if(A.inCoeffDomain())
293 {
294 tryInvert( A, M, P, fail );
295 if(fail)
296 return;
297 result = 1;
298 return;
299 }
300 if(B.inCoeffDomain())
301 {
302 tryInvert( B, M, P, fail );
303 if(fail)
304 return;
305 result = 1;
306 return;
307 }
308 // here: both not inCoeffDomain
309 if( A.degree() > B.degree() )
310 {
311 P = A; result = B;
312 }
313 else
314 {
315 P = B; result = A;
316 }
317 CanonicalForm inv;
318 if( result.isZero() )
319 {
320 tryInvert( Lc(P), M, inv, fail );
321 if(fail)
322 return;
323 result = inv*P; // monify result (not reduced, yet)
324 result= reduce (result, M);
325 return;
326 }
327 Variable x = P.mvar();
329 // here: degree(P) >= degree(result)
330 while(true)
331 {
332 tryDivrem (P, result, Q, rem, inv, M, fail);
333 if (fail)
334 return;
335 if( rem.isZero() )
336 {
337 result *= inv;
338 result= reduce (result, M);
339 return;
340 }
341 if(result.degree(x) >= rem.degree(x))
342 {
343 P = result;
344 result = rem;
345 }
346 else
347 P = rem;
348 }
349}
350#endif
351
352static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
353static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
354static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail );
355
356static inline CanonicalForm
358 const CanonicalForm & newtonPoly, const CanonicalForm & oldInterPoly,
359 const Variable & x, const CanonicalForm& M, bool& fail)
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}
371
372static void leadDeg(const CanonicalForm & f, int degs[])
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}
385
386void tryBrownGCD( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, bool & fail, bool topLevel )
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;
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}
669
670static CanonicalForm
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}
720
722{
723#if defined(HAVE_NTL) || defined(HAVE_FLINT)
724 return myicontent( f, 0 );
725#else
726 return 1;
727#endif
728}
729
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;
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 {
877 newtonDivrem (f, tmp, Q, R);
878 if (R.isZero())
879 {
880 newtonDivrem (g, tmp, Q, R);
881 if (R.isZero())
882 {
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}
932
933
934
935bool isLess(int *a, int *b, int lower, int upper)
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}
944
945
946bool isEqual(int *a, int *b, int lower, int upper)
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}
953
954
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}
962
963#ifndef HAVE_NTL
964void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail )
965{ // F, G are univariate polynomials (i.e. they have exactly one polynomial variable)
966 // F and G must have the same level AND level > 0
967 // we try to calculate gcd(F,G) = s*F + t*G
968 // if a zero divisor is encountered, 'fail' is set to one
969 // M is assumed to be monic
971 if(F.inCoeffDomain())
972 {
973 tryInvert( F, M, P, fail );
974 if(fail)
975 return;
976 result = 1;
977 s = P; t = 0;
978 return;
979 }
980 if(G.inCoeffDomain())
981 {
982 tryInvert( G, M, P, fail );
983 if(fail)
984 return;
985 result = 1;
986 s = 0; t = P;
987 return;
988 }
989 // here: both not inCoeffDomain
990 CanonicalForm inv, rem, tmp, u, v, q, sum=0;
991 if( F.degree() > G.degree() )
992 {
993 P = F; result = G; s=v=0; t=u=1;
994 }
995 else
996 {
997 P = G; result = F; s=v=1; t=u=0;
998 }
999 Variable x = P.mvar();
1000 // here: degree(P) >= degree(result)
1001 while(true)
1002 {
1003 tryDivrem (P, result, q, rem, inv, M, fail);
1004 if(fail)
1005 return;
1006 if( rem.isZero() )
1007 {
1008 s*=inv;
1009 s= reduce (s, M);
1010 t*=inv;
1011 t= reduce (t, M);
1012 result *= inv; // monify result
1013 result= reduce (result, M);
1014 return;
1015 }
1016 sum += q;
1017 if(result.degree(x) >= rem.degree(x))
1018 {
1019 P=result;
1020 result=rem;
1021 tmp=u-sum*s;
1022 u=s;
1023 s=tmp;
1024 tmp=v-sum*t;
1025 v=t;
1026 t=tmp;
1027 sum = 0; // reset
1028 }
1029 else
1030 P = rem;
1031 }
1032}
1033#endif
1034
1035static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
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}
1045
1046
1047static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
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}
1064
1065
1066static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail )
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}
1082
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
This file defines functions for conversion to FLINT (www.flintlib.org) and back.
Rational abs(const Rational &a)
Definition: GMPrat.cc:436
ZZX convertFacCF2NTLZZX(const CanonicalForm &f)
Definition: NTLconvert.cc:691
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
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 )
bool isOn(int sw)
switches
void On(int sw)
switches
void Off(int sw)
switches
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Header for factory's main class CanonicalForm.
CanonicalForm mapinto(const CanonicalForm &f)
CanonicalForm lc(const CanonicalForm &f)
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
int * degrees(const CanonicalForm &f, int *degs=0)
int * degrees ( const CanonicalForm & f, int * degs )
Definition: cf_ops.cc:493
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
int degree(const CanonicalForm &f)
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 FACTORY_PUBLIC swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition: cf_ops.cc:168
CanonicalForm Lc(const CanonicalForm &f)
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
CanonicalForm LC(const CanonicalForm &f)
int FACTORY_PUBLIC getCharacteristic()
Definition: cf_char.cc:70
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
const CanonicalForm CFMap CFMap & N
Definition: cfGcdAlgExt.cc:56
bool isLess(int *a, int *b, int lower, int upper)
Definition: cfGcdAlgExt.cc:935
int both_non_zero
Definition: cfGcdAlgExt.cc:68
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
static void leadDeg(const CanonicalForm &f, int degs[])
Definition: cfGcdAlgExt.cc:372
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
int * degsf
Definition: cfGcdAlgExt.cc:59
static CanonicalForm trycontent(const CanonicalForm &f, const Variable &x, const CanonicalForm &M, bool &fail)
int f_zero
Definition: cfGcdAlgExt.cc:69
const CanonicalForm CFMap CFMap bool topLevel
Definition: cfGcdAlgExt.cc:57
bool isEqual(int *a, int *b, int lower, int upper)
Definition: cfGcdAlgExt.cc:946
static CanonicalForm myicontent(const CanonicalForm &f, const CanonicalForm &c)
Definition: cfGcdAlgExt.cc:671
int both_zero
Definition: cfGcdAlgExt.cc:71
int Flevel
Definition: cfGcdAlgExt.cc:72
static CanonicalForm tryvcontent(const CanonicalForm &f, const Variable &x, const CanonicalForm &M, bool &fail)
const CanonicalForm & G
Definition: cfGcdAlgExt.cc:55
int Glevel
Definition: cfGcdAlgExt.cc:73
void tryInvert(const CanonicalForm &F, const CanonicalForm &M, CanonicalForm &inv, bool &fail)
Definition: cfGcdAlgExt.cc:221
static CanonicalForm trycf_content(const CanonicalForm &f, const CanonicalForm &g, const CanonicalForm &M, bool &fail)
const CanonicalForm CFMap & M
Definition: cfGcdAlgExt.cc:55
int g_zero
Definition: cfGcdAlgExt.cc:70
int * degsg
Definition: cfGcdAlgExt.cc:60
CanonicalForm QGCD(const CanonicalForm &F, const CanonicalForm &G)
gcd over Q(a)
Definition: cfGcdAlgExt.cc:730
CanonicalForm firstLC(const CanonicalForm &f)
Definition: cfGcdAlgExt.cc:955
GCD over Q(a)
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
int p
Definition: cfModGcd.cc:4078
cl
Definition: cfModGcd.cc:4100
g
Definition: cfModGcd.cc:4090
CanonicalForm gcdcfcg
Definition: cfModGcd.cc:4101
CanonicalForm cf
Definition: cfModGcd.cc:4083
bool equal
Definition: cfModGcd.cc:4126
CanonicalForm test
Definition: cfModGcd.cc:4096
CanonicalForm b
Definition: cfModGcd.cc:4103
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
This file defines functions for univariate GCD and extended GCD over Z/p[t]/(f)[x] for reducible f.
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
univariate Gcd over finite fields and Z, extended GCD over finite fields and Q
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 )
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
declarations of higher level algorithms.
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
assertions for Factory
#define ASSERT(expression, message)
Definition: cf_assert.h:99
factory switches.
static const int SW_USE_QGCD
set to 1 to use Encarnacion GCD over Q(a)
Definition: cf_defs.h:43
#define DELETE_ARRAY(P)
Definition: cf_defs.h:65
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:31
#define NEW_ARRAY(T, N)
Definition: cf_defs.h:64
generate integers, elements of finite fields
Iterators for CanonicalForm's.
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
map polynomials
int cf_getBigPrime(int i)
Definition: cf_primes.cc:39
int cf_getNumBigPrimes()
Definition: cf_primes.cc:45
access to prime tables
FILE * f
Definition: checklibs.c:9
class to iterate through CanonicalForm's
Definition: cf_iter.h:44
class CFMap
Definition: cf_map.h:85
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.
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 inCoeffDomain() const
int level() const
level() returns the level of CO.
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
bool inBaseDomain() const
CanonicalForm lc() const
CanonicalForm CanonicalForm::lc (), Lc (), LC (), LC ( v ) const.
bool isUnivariate() const
generate all elements in F_p starting from 0
Definition: cf_generator.h:56
bool hasItems() const
Definition: cf_generator.cc:35
CanonicalForm item() const
Definition: cf_generator.cc:40
factory's class for variables
Definition: variable.h:33
int level() const
Definition: variable.h:49
Variable alpha
Definition: facAbsBiFact.cc:51
return result
Definition: facAbsBiFact.cc:75
const CanonicalForm int s
Definition: facAbsFact.cc:51
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53
CanonicalForm mipo
Definition: facAlgExt.cc:57
CanonicalForm alg_content(const CanonicalForm &f, const CFList &as)
Definition: facAlgFunc.cc:42
b *CanonicalForm B
Definition: facBivar.cc:52
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
int j
Definition: facHensel.cc:110
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
This file defines functions for fast multiplication and division with remainder.
bool isZero(const CFArray &A)
checks if entries of A are zero
#define const
Definition: fegetopt.c:39
static number Farey(number, number, const coeffs)
Definition: flintcf_Q.cc:445
some useful template functions.
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
#define D(A)
Definition: gentable.cc:131
void rem(unsigned long *a, unsigned long *q, unsigned long p, int &dega, int degq)
Definition: minpoly.cc:572
#define R
Definition: sirandom.c:27
#define A
Definition: sirandom.c:24
#define Q
Definition: sirandom.c:26
#define TIMING_DEFINE_PRINT(t)
Definition: timing.h:95
#define TIMING_START(t)
Definition: timing.h:92
#define TIMING_END_AND_PRINT(t, msg)
Definition: timing.h:94
void setReduce(const Variable &alpha, bool reduce)
Definition: variable.cc:238
bool getReduce(const Variable &alpha)
Definition: variable.cc:232
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition: variable.cc:207
int gcd(int a, int b)
Definition: walkSupport.cc:836