My Project
Loading...
Searching...
No Matches
cfModGcd.cc
Go to the documentation of this file.
1// -*- c++ -*-
2//*****************************************************************************
3/** @file cfModGcd.cc
4 *
5 * This file implements the GCD of two polynomials over \f$ F_{p} \f$ ,
6 * \f$ F_{p}(\alpha ) \f$, GF or Z based on Alg. 7.1. and 7.2. as described in
7 * "Algorithms for Computer Algebra" by Geddes, Czapor, Labahn via modular
8 * computations. And sparse modular variants as described in Zippel
9 * "Effective Polynomial Computation", deKleine, Monagan, Wittkopf
10 * "Algorithms for the non-monic case of the sparse modular GCD algorithm" and
11 * Javadi "A new solution to the normalization problem"
12 *
13 * @author Martin Lee
14 * @date 22.10.2009
15 *
16 * @par Copyright:
17 * (c) by The SINGULAR Team, see LICENSE file
18 *
19**/
20//*****************************************************************************
21
22
23#include "config.h"
24
25
26#include "cf_assert.h"
27#include "debug.h"
28#include "timing.h"
29
30#include "canonicalform.h"
31#include "cfGcdUtil.h"
32#include "cf_map.h"
33#include "cf_util.h"
34#include "cf_irred.h"
36#include "cf_random.h"
37#include "cf_reval.h"
38#include "facHensel.h"
39#include "cf_iter.h"
40#include "cfNewtonPolygon.h"
41#include "cf_algorithm.h"
42#include "cf_primes.h"
43
44// inline helper functions:
45#include "cf_map_ext.h"
46
47#ifdef HAVE_NTL
48#include "NTLconvert.h"
49#endif
50
51#ifdef HAVE_FLINT
52#include "FLINTconvert.h"
53#endif
54
55#include "cfModGcd.h"
56
57TIMING_DEFINE_PRINT(gcd_recursion)
58TIMING_DEFINE_PRINT(newton_interpolation)
59TIMING_DEFINE_PRINT(termination_test)
60TIMING_DEFINE_PRINT(ez_p_compress)
61TIMING_DEFINE_PRINT(ez_p_hensel_lift)
62TIMING_DEFINE_PRINT(ez_p_eval)
63TIMING_DEFINE_PRINT(ez_p_content)
64
65bool
69{
70 CanonicalForm LCCand= abs (LC (cand));
71 if (LCCand*abs (LC (coF)) == abs (LC (F)))
72 {
73 if (LCCand*abs (LC (coG)) == abs (LC (G)))
74 {
75 if (abs (cand)*abs (coF) == abs (F))
76 {
77 if (abs (cand)*abs (coG) == abs (G))
78 return true;
79 }
80 return false;
81 }
82 return false;
83 }
84 return false;
85}
86
87#if defined(HAVE_NTL)|| defined(HAVE_FLINT)
88
89/// compressing two polynomials F and G, M is used for compressing,
90/// N to reverse the compression
91int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
92 CFMap & N, bool topLevel)
93{
94 int n= tmax (F.level(), G.level());
95 int * degsf= NEW_ARRAY(int,n + 1);
96 int * degsg= NEW_ARRAY(int,n + 1);
97
98 for (int i = n; i >= 0; i--)
99 degsf[i]= degsg[i]= 0;
100
101 degsf= degrees (F, degsf);
102 degsg= degrees (G, degsg);
103
104 int both_non_zero= 0;
105 int f_zero= 0;
106 int g_zero= 0;
107 int both_zero= 0;
108
109 if (topLevel)
110 {
111 for (int i= 1; i <= n; i++)
112 {
113 if (degsf[i] != 0 && degsg[i] != 0)
114 {
116 continue;
117 }
118 if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
119 {
120 f_zero++;
121 continue;
122 }
123 if (degsg[i] == 0 && degsf[i] && i <= F.level())
124 {
125 g_zero++;
126 continue;
127 }
128 }
129
130 if (both_non_zero == 0)
131 {
134 return 0;
135 }
136
137 // map Variables which do not occur in both polynomials to higher levels
138 int k= 1;
139 int l= 1;
140 for (int i= 1; i <= n; i++)
141 {
142 if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
143 {
144 if (k + both_non_zero != i)
145 {
146 M.newpair (Variable (i), Variable (k + both_non_zero));
147 N.newpair (Variable (k + both_non_zero), Variable (i));
148 }
149 k++;
150 }
151 if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
152 {
153 if (l + g_zero + both_non_zero != i)
154 {
155 M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
156 N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
157 }
158 l++;
159 }
160 }
161
162 // sort Variables x_{i} in increasing order of
163 // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
164 int m= tmax (F.level(), G.level());
165 int min_max_deg;
167 l= 0;
168 int i= 1;
169 while (k > 0)
170 {
171 if (degsf [i] != 0 && degsg [i] != 0)
172 min_max_deg= tmax (degsf[i], degsg[i]);
173 else
174 min_max_deg= 0;
175 while (min_max_deg == 0)
176 {
177 i++;
178 if (degsf [i] != 0 && degsg [i] != 0)
179 min_max_deg= tmax (degsf[i], degsg[i]);
180 else
181 min_max_deg= 0;
182 }
183 for (int j= i + 1; j <= m; j++)
184 {
185 if (degsf[j] != 0 && degsg [j] != 0 &&
186 tmax (degsf[j],degsg[j]) <= min_max_deg)
187 {
188 min_max_deg= tmax (degsf[j], degsg[j]);
189 l= j;
190 }
191 }
192 if (l != 0)
193 {
194 if (l != k)
195 {
196 M.newpair (Variable (l), Variable(k));
197 N.newpair (Variable (k), Variable(l));
198 degsf[l]= 0;
199 degsg[l]= 0;
200 l= 0;
201 }
202 else
203 {
204 degsf[l]= 0;
205 degsg[l]= 0;
206 l= 0;
207 }
208 }
209 else if (l == 0)
210 {
211 if (i != k)
212 {
213 M.newpair (Variable (i), Variable (k));
214 N.newpair (Variable (k), Variable (i));
215 degsf[i]= 0;
216 degsg[i]= 0;
217 }
218 else
219 {
220 degsf[i]= 0;
221 degsg[i]= 0;
222 }
223 i++;
224 }
225 k--;
226 }
227 }
228 else
229 {
230 //arrange Variables such that no gaps occur
231 for (int i= 1; i <= n; i++)
232 {
233 if (degsf[i] == 0 && degsg[i] == 0)
234 {
235 both_zero++;
236 continue;
237 }
238 else
239 {
240 if (both_zero != 0)
241 {
242 M.newpair (Variable (i), Variable (i - both_zero));
243 N.newpair (Variable (i - both_zero), Variable (i));
244 }
245 }
246 }
247 }
248
251
252 return 1;
253}
254
255static inline CanonicalForm
256uni_content (const CanonicalForm & F);
257
260{
261 if (F.inCoeffDomain())
262 return F.genOne();
263 if (F.level() == x.level() && F.isUnivariate())
264 return F;
265 if (F.level() != x.level() && F.isUnivariate())
266 return F.genOne();
267
268 if (x.level() != 1)
269 {
270 CanonicalForm f= swapvar (F, x, Variable (1));
272 return swapvar (result, x, Variable (1));
273 }
274 else
275 return uni_content (F);
276}
277
278/// compute the content of F, where F is considered as an element of
279/// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$
280static inline CanonicalForm
282{
283 if (F.inBaseDomain())
284 return F.genOne();
285 if (F.level() == 1 && F.isUnivariate())
286 return F;
287 if (F.level() != 1 && F.isUnivariate())
288 return F.genOne();
289 if (degree (F,1) == 0) return F.genOne();
290
291 int l= F.level();
292 if (l == 2)
293 return content(F);
294 else
295 {
296 CanonicalForm pol, c = 0;
297 CFIterator i = F;
298 for (; i.hasTerms(); i++)
299 {
300 pol= i.coeff();
301 pol= uni_content (pol);
302 c= gcd (c, pol);
303 if (c.isOne())
304 return c;
305 }
306 return c;
307 }
308}
309
312 CanonicalForm& contentF, CanonicalForm& contentG,
313 CanonicalForm& ppF, CanonicalForm& ppG, const int d)
314{
315 CanonicalForm uniContentF, uniContentG, gcdcFcG;
316 contentF= 1;
317 contentG= 1;
318 ppF= F;
319 ppG= G;
321 for (int i= 1; i <= d; i++)
322 {
323 uniContentF= uni_content (F, Variable (i));
324 uniContentG= uni_content (G, Variable (i));
325 gcdcFcG= gcd (uniContentF, uniContentG);
326 contentF *= uniContentF;
327 contentG *= uniContentG;
328 ppF /= uniContentF;
329 ppG /= uniContentG;
330 result *= gcdcFcG;
331 }
332 return result;
333}
334
335/// compute the leading coefficient of F, where F is considered as an element
336/// of \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$, order on
337/// \f$ Mon (x_{2},\ldots ,x_{n}) \f$ is dp.
338static inline
340{
341 if (F.level() > 1)
342 {
343 Variable x= Variable (2);
344 int deg= totaldegree (F, x, F.mvar());
345 for (CFIterator i= F; i.hasTerms(); i++)
346 {
347 if (i.exp() + totaldegree (i.coeff(), x, i.coeff().mvar()) == deg)
348 return uni_lcoeff (i.coeff());
349 }
350 }
351 return F;
352}
353
354/// Newton interpolation - Incremental algorithm.
355/// Given a list of values alpha_i and a list of polynomials u_i, 1 <= i <= n,
356/// computes the interpolation polynomial assuming that
357/// the polynomials in u are the results of evaluating the variabe x
358/// of the unknown polynomial at the alpha values.
359/// This incremental version receives only the values of alpha_n and u_n and
360/// the previous interpolation polynomial for points 1 <= i <= n-1, and computes
361/// the polynomial interpolating in all the points.
362/// newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1})
363static inline CanonicalForm
365 const CanonicalForm & newtonPoly, const CanonicalForm & oldInterPoly,
366 const Variable & x)
367{
368 CanonicalForm interPoly;
369
370 interPoly= oldInterPoly + ((u - oldInterPoly(alpha, x))/newtonPoly(alpha, x))
371 *newtonPoly;
372 return interPoly;
373}
374
375/// compute a random element a of \f$ F_{p}(\alpha ) \f$ ,
376/// s.t. F(a) \f$ \neq 0 \f$ , F is a univariate polynomial, returns
377/// fail if there are no field elements left which have not been used before
378static inline CanonicalForm
379randomElement (const CanonicalForm & F, const Variable & alpha, CFList & list,
380 bool & fail)
381{
382 fail= false;
383 Variable x= F.mvar();
384 AlgExtRandomF genAlgExt (alpha);
385 FFRandom genFF;
386 CanonicalForm random, mipo;
387 mipo= getMipo (alpha);
388 int p= getCharacteristic ();
389 int d= degree (mipo);
390 double bound= pow ((double) p, (double) d);
391 do
392 {
393 if (list.length() == bound)
394 {
395 fail= true;
396 break;
397 }
398 if (list.length() < p)
399 {
400 random= genFF.generate();
401 while (find (list, random))
402 random= genFF.generate();
403 }
404 else
405 {
406 random= genAlgExt.generate();
407 while (find (list, random))
408 random= genAlgExt.generate();
409 }
410 if (F (random, x) == 0)
411 {
412 list.append (random);
413 continue;
414 }
415 } while (find (list, random));
416 return random;
417}
418
419static inline
421{
422 int i, m;
423 // extension of F_p needed
424 if (alpha.level() == 1)
425 {
426 i= 1;
427 m= 2;
428 } //extension of F_p(alpha)
429 if (alpha.level() != 1)
430 {
431 i= 4;
432 m= degree (getMipo (alpha));
433 }
434 #ifdef HAVE_FLINT
435 nmod_poly_t Irredpoly;
437 nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, i*m+1);
438 CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
439 nmod_poly_clear(Irredpoly);
440 #else
442 {
444 zz_p::init (getCharacteristic());
445 }
446 zz_pX NTLIrredpoly;
447 BuildIrred (NTLIrredpoly, i*m);
448 CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
449 #endif
450 return rootOf (newMipo);
451}
452
453#if defined(HAVE_NTL) || defined(HAVE_FLINT)
455modGCDFq (const CanonicalForm& F, const CanonicalForm& G,
457 Variable & alpha, CFList& l, bool& topLevel);
458#endif
459
460#if defined(HAVE_NTL) || defined(HAVE_FLINT)
463 Variable & alpha, CFList& l, bool& topLevel)
464{
465 CanonicalForm dummy1, dummy2;
466 CanonicalForm result= modGCDFq (F, G, dummy1, dummy2, alpha, l,
467 topLevel);
468 return result;
469}
470#endif
471
472/// GCD of F and G over \f$ F_{p}(\alpha ) \f$ ,
473/// l and topLevel are only used internally, output is monic
474/// based on Alg. 7.2. as described in "Algorithms for
475/// Computer Algebra" by Geddes, Czapor, Labahn
476#if defined(HAVE_NTL) || defined(HAVE_FLINT)
480 Variable & alpha, CFList& l, bool& topLevel)
481{
482 CanonicalForm A= F;
484 if (F.isZero() && degree(G) > 0)
485 {
486 coF= 0;
487 coG= Lc (G);
488 return G/Lc(G);
489 }
490 else if (G.isZero() && degree (F) > 0)
491 {
492 coF= Lc (F);
493 coG= 0;
494 return F/Lc(F);
495 }
496 else if (F.isZero() && G.isZero())
497 {
498 coF= coG= 0;
499 return F.genOne();
500 }
501 if (F.inBaseDomain() || G.inBaseDomain())
502 {
503 coF= F;
504 coG= G;
505 return F.genOne();
506 }
507 if (F.isUnivariate() && fdivides(F, G, coG))
508 {
509 coF= Lc (F);
510 return F/Lc(F);
511 }
512 if (G.isUnivariate() && fdivides(G, F, coF))
513 {
514 coG= Lc (G);
515 return G/Lc(G);
516 }
517 if (F == G)
518 {
519 coF= coG= Lc (F);
520 return F/Lc(F);
521 }
522
523 CFMap M,N;
524 int best_level= myCompress (A, B, M, N, topLevel);
525
526 if (best_level == 0)
527 {
528 coF= F;
529 coG= G;
530 return B.genOne();
531 }
532
533 A= M(A);
534 B= M(B);
535
536 Variable x= Variable(1);
537
538 //univariate case
539 if (A.isUnivariate() && B.isUnivariate())
540 {
542 coF= N (A/result);
543 coG= N (B/result);
544 return N (result);
545 }
546
547 CanonicalForm cA, cB; // content of A and B
548 CanonicalForm ppA, ppB; // primitive part of A and B
549 CanonicalForm gcdcAcB;
550
551 cA = uni_content (A);
552 cB = uni_content (B);
553 gcdcAcB= gcd (cA, cB);
554 ppA= A/cA;
555 ppB= B/cB;
556
557 CanonicalForm lcA, lcB; // leading coefficients of A and B
558 CanonicalForm gcdlcAlcB;
559
560 lcA= uni_lcoeff (ppA);
561 lcB= uni_lcoeff (ppB);
562
563 gcdlcAlcB= gcd (lcA, lcB);
564
565 int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
566
567 if (d == 0)
568 {
569 coF= N (ppA*(cA/gcdcAcB));
570 coG= N (ppB*(cB/gcdcAcB));
571 return N(gcdcAcB);
572 }
573
574 int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
575 if (d0 < d)
576 d= d0;
577 if (d == 0)
578 {
579 coF= N (ppA*(cA/gcdcAcB));
580 coG= N (ppB*(cB/gcdcAcB));
581 return N(gcdcAcB);
582 }
583
584 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
585 CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
586 coG_m, ppCoF, ppCoG;
587
588 newtonPoly= 1;
589 m= gcdlcAlcB;
590 G_m= 0;
591 coF= 0;
592 coG= 0;
593 H= 0;
594 bool fail= false;
595 topLevel= false;
596 bool inextension= false;
597 Variable V_buf= alpha, V_buf4= alpha;
598 CanonicalForm prim_elem, im_prim_elem;
599 CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
600 CFList source, dest;
601 int bound1= degree (ppA, 1);
602 int bound2= degree (ppB, 1);
603 do
604 {
605 random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
606 if (fail)
607 {
608 source= CFList();
609 dest= CFList();
610
611 Variable V_buf3= V_buf;
612 V_buf= chooseExtension (V_buf);
613 bool prim_fail= false;
614 Variable V_buf2;
615 prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
616 if (V_buf4 == alpha)
617 prim_elem_alpha= prim_elem;
618
619 if (V_buf3 != V_buf4)
620 {
621 m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
622 G_m= mapDown (G_m, prim_elem, im_prim_elem, V_buf4, source, dest);
623 coF_m= mapDown (coF_m, prim_elem, im_prim_elem, V_buf4, source, dest);
624 coG_m= mapDown (coG_m, prim_elem, im_prim_elem, V_buf4, source, dest);
625 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
626 source, dest);
627 ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
628 ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
629 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
630 source, dest);
631 lcA= mapDown (lcA, prim_elem, im_prim_elem, V_buf4, source, dest);
632 lcB= mapDown (lcB, prim_elem, im_prim_elem, V_buf4, source, dest);
633 for (CFListIterator i= l; i.hasItem(); i++)
634 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
635 source, dest);
636 }
637
638 ASSERT (!prim_fail, "failure in integer factorizer");
639 if (prim_fail)
640 ; //ERROR
641 else
642 im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
643
644 if (V_buf4 == alpha)
645 im_prim_elem_alpha= im_prim_elem;
646 else
647 im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
648 im_prim_elem, source, dest);
649 DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
650 DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
651 inextension= true;
652 for (CFListIterator i= l; i.hasItem(); i++)
653 i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
654 im_prim_elem, source, dest);
655 m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
656 G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
657 coF_m= mapUp (coF_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
658 coG_m= mapUp (coG_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
659 newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
660 source, dest);
661 ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
662 ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
663 gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
664 source, dest);
665 lcA= mapUp (lcA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
666 lcB= mapUp (lcB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
667
668 fail= false;
669 random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
670 DEBOUTLN (cerr, "fail= " << fail);
671 CFList list;
672 TIMING_START (gcd_recursion);
673 G_random_element=
674 modGCDFq (ppA (random_element, x), ppB (random_element, x),
675 coF_random_element, coG_random_element, V_buf,
676 list, topLevel);
677 TIMING_END_AND_PRINT (gcd_recursion,
678 "time for recursive call: ");
679 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
680 V_buf4= V_buf;
681 }
682 else
683 {
684 CFList list;
685 TIMING_START (gcd_recursion);
686 G_random_element=
687 modGCDFq (ppA(random_element, x), ppB(random_element, x),
688 coF_random_element, coG_random_element, V_buf,
689 list, topLevel);
690 TIMING_END_AND_PRINT (gcd_recursion,
691 "time for recursive call: ");
692 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
693 }
694
695 if (!G_random_element.inCoeffDomain())
696 d0= totaldegree (G_random_element, Variable(2),
697 Variable (G_random_element.level()));
698 else
699 d0= 0;
700
701 if (d0 == 0)
702 {
703 if (inextension)
704 {
705 CFList u, v;
706 ppA= mapDown (ppA, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
707 ppB= mapDown (ppB, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
708 prune1 (alpha);
709 }
710 coF= N (ppA*(cA/gcdcAcB));
711 coG= N (ppB*(cB/gcdcAcB));
712 return N(gcdcAcB);
713 }
714 if (d0 > d)
715 {
716 if (!find (l, random_element))
717 l.append (random_element);
718 continue;
719 }
720
721 G_random_element=
722 (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
723 * G_random_element;
724 coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
725 *coF_random_element;
726 coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
727 *coG_random_element;
728
729 if (!G_random_element.inCoeffDomain())
730 d0= totaldegree (G_random_element, Variable(2),
731 Variable (G_random_element.level()));
732 else
733 d0= 0;
734
735 if (d0 < d)
736 {
737 m= gcdlcAlcB;
738 newtonPoly= 1;
739 G_m= 0;
740 d= d0;
741 coF_m= 0;
742 coG_m= 0;
743 }
744
745 TIMING_START (newton_interpolation);
746 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
747 coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
748 coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
749 TIMING_END_AND_PRINT (newton_interpolation,
750 "time for newton interpolation: ");
751
752 //termination test
753 if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
754 {
755 TIMING_START (termination_test);
756 if (gcdlcAlcB.isOne())
757 cH= 1;
758 else
759 cH= uni_content (H);
760 ppH= H/cH;
761 ppH /= Lc (ppH);
762 CanonicalForm lcppH= gcdlcAlcB/cH;
763 CanonicalForm ccoF= lcppH/Lc (lcppH);
764 CanonicalForm ccoG= lcppH/Lc (lcppH);
765 ppCoF= coF/ccoF;
766 ppCoG= coG/ccoG;
767 if (inextension)
768 {
769 if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
770 (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
771 terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
772 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
773 {
774 CFList u, v;
775 DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
776 ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
777 ppCoF= mapDown (ppCoF, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
778 ppCoG= mapDown (ppCoG, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
779 DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
780 coF= N ((cA/gcdcAcB)*ppCoF);
781 coG= N ((cB/gcdcAcB)*ppCoG);
782 TIMING_END_AND_PRINT (termination_test,
783 "time for successful termination test Fq: ");
784 prune1 (alpha);
785 return N(gcdcAcB*ppH);
786 }
787 }
788 else if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
789 (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
790 terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
791 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
792 {
793 coF= N ((cA/gcdcAcB)*ppCoF);
794 coG= N ((cB/gcdcAcB)*ppCoG);
795 TIMING_END_AND_PRINT (termination_test,
796 "time for successful termination test Fq: ");
797 return N(gcdcAcB*ppH);
798 }
799 TIMING_END_AND_PRINT (termination_test,
800 "time for unsuccessful termination test Fq: ");
801 }
802
803 G_m= H;
804 coF_m= coF;
805 coG_m= coG;
806 newtonPoly= newtonPoly*(x - random_element);
807 m= m*(x - random_element);
808 if (!find (l, random_element))
809 l.append (random_element);
810 } while (1);
811}
812#endif
813
814/// compute a random element a of GF, s.t. F(a) \f$ \neq 0 \f$ , F is a
815/// univariate polynomial, returns fail if there are no field elements left
816/// which have not been used before
817static inline
819GFRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
820{
821 fail= false;
822 Variable x= F.mvar();
823 GFRandom genGF;
824 CanonicalForm random;
825 int p= getCharacteristic();
826 int d= getGFDegree();
827 int bound= ipower (p, d);
828 do
829 {
830 if (list.length() == bound)
831 {
832 fail= true;
833 break;
834 }
835 if (list.length() < 1)
836 random= 0;
837 else
838 {
839 random= genGF.generate();
840 while (find (list, random))
841 random= genGF.generate();
842 }
843 if (F (random, x) == 0)
844 {
845 list.append (random);
846 continue;
847 }
848 } while (find (list, random));
849 return random;
850}
851
853modGCDGF (const CanonicalForm& F, const CanonicalForm& G,
855 CFList& l, bool& topLevel);
856
859 bool& topLevel)
860{
861 CanonicalForm dummy1, dummy2;
862 CanonicalForm result= modGCDGF (F, G, dummy1, dummy2, l, topLevel);
863 return result;
864}
865
866/// GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for
867/// Computer Algebra" by Geddes, Czapor, Labahn
868/// Usually this algorithm will be faster than modGCDFq since GF has
869/// faster field arithmetics, however it might fail if the input is large since
870/// the size of the base field is bounded by 2^16, output is monic
874 CFList& l, bool& topLevel)
875{
876 CanonicalForm A= F;
878 if (F.isZero() && degree(G) > 0)
879 {
880 coF= 0;
881 coG= Lc (G);
882 return G/Lc(G);
883 }
884 else if (G.isZero() && degree (F) > 0)
885 {
886 coF= Lc (F);
887 coG= 0;
888 return F/Lc(F);
889 }
890 else if (F.isZero() && G.isZero())
891 {
892 coF= coG= 0;
893 return F.genOne();
894 }
895 if (F.inBaseDomain() || G.inBaseDomain())
896 {
897 coF= F;
898 coG= G;
899 return F.genOne();
900 }
901 if (F.isUnivariate() && fdivides(F, G, coG))
902 {
903 coF= Lc (F);
904 return F/Lc(F);
905 }
906 if (G.isUnivariate() && fdivides(G, F, coF))
907 {
908 coG= Lc (G);
909 return G/Lc(G);
910 }
911 if (F == G)
912 {
913 coF= coG= Lc (F);
914 return F/Lc(F);
915 }
916
917 CFMap M,N;
918 int best_level= myCompress (A, B, M, N, topLevel);
919
920 if (best_level == 0)
921 {
922 coF= F;
923 coG= G;
924 return B.genOne();
925 }
926
927 A= M(A);
928 B= M(B);
929
930 Variable x= Variable(1);
931
932 //univariate case
933 if (A.isUnivariate() && B.isUnivariate())
934 {
936 coF= N (A/result);
937 coG= N (B/result);
938 return N (result);
939 }
940
941 CanonicalForm cA, cB; // content of A and B
942 CanonicalForm ppA, ppB; // primitive part of A and B
943 CanonicalForm gcdcAcB;
944
945 cA = uni_content (A);
946 cB = uni_content (B);
947 gcdcAcB= gcd (cA, cB);
948 ppA= A/cA;
949 ppB= B/cB;
950
951 CanonicalForm lcA, lcB; // leading coefficients of A and B
952 CanonicalForm gcdlcAlcB;
953
954 lcA= uni_lcoeff (ppA);
955 lcB= uni_lcoeff (ppB);
956
957 gcdlcAlcB= gcd (lcA, lcB);
958
959 int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
960 if (d == 0)
961 {
962 coF= N (ppA*(cA/gcdcAcB));
963 coG= N (ppB*(cB/gcdcAcB));
964 return N(gcdcAcB);
965 }
966 int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
967 if (d0 < d)
968 d= d0;
969 if (d == 0)
970 {
971 coF= N (ppA*(cA/gcdcAcB));
972 coG= N (ppB*(cB/gcdcAcB));
973 return N(gcdcAcB);
974 }
975
976 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
977 CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
978 coG_m, ppCoF, ppCoG;
979
980 newtonPoly= 1;
981 m= gcdlcAlcB;
982 G_m= 0;
983 coF= 0;
984 coG= 0;
985 H= 0;
986 bool fail= false;
987 topLevel= false;
988 bool inextension= false;
989 int p=-1;
990 int k= getGFDegree();
991 int kk;
992 int expon;
993 char gf_name_buf= gf_name;
994 int bound1= degree (ppA, 1);
995 int bound2= degree (ppB, 1);
996 do
997 {
998 random_element= GFRandomElement (m*lcA*lcB, l, fail);
999 if (fail)
1000 {
1002 expon= 2;
1003 kk= getGFDegree();
1004 if (ipower (p, kk*expon) < (1 << 16))
1005 setCharacteristic (p, kk*(int)expon, 'b');
1006 else
1007 {
1008 expon= (int) floor((log ((double)((1<<16) - 1)))/(log((double)p)*kk));
1009 ASSERT (expon >= 2, "not enough points in modGCDGF");
1010 setCharacteristic (p, (int)(kk*expon), 'b');
1011 }
1012 inextension= true;
1013 fail= false;
1014 for (CFListIterator i= l; i.hasItem(); i++)
1015 i.getItem()= GFMapUp (i.getItem(), kk);
1016 m= GFMapUp (m, kk);
1017 G_m= GFMapUp (G_m, kk);
1018 newtonPoly= GFMapUp (newtonPoly, kk);
1019 coF_m= GFMapUp (coF_m, kk);
1020 coG_m= GFMapUp (coG_m, kk);
1021 ppA= GFMapUp (ppA, kk);
1022 ppB= GFMapUp (ppB, kk);
1023 gcdlcAlcB= GFMapUp (gcdlcAlcB, kk);
1024 lcA= GFMapUp (lcA, kk);
1025 lcB= GFMapUp (lcB, kk);
1026 random_element= GFRandomElement (m*lcA*lcB, l, fail);
1027 DEBOUTLN (cerr, "fail= " << fail);
1028 CFList list;
1029 TIMING_START (gcd_recursion);
1030 G_random_element= modGCDGF (ppA(random_element, x), ppB(random_element, x),
1031 coF_random_element, coG_random_element,
1032 list, topLevel);
1033 TIMING_END_AND_PRINT (gcd_recursion,
1034 "time for recursive call: ");
1035 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1036 }
1037 else
1038 {
1039 CFList list;
1040 TIMING_START (gcd_recursion);
1041 G_random_element= modGCDGF (ppA(random_element, x), ppB(random_element, x),
1042 coF_random_element, coG_random_element,
1043 list, topLevel);
1044 TIMING_END_AND_PRINT (gcd_recursion,
1045 "time for recursive call: ");
1046 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1047 }
1048
1049 if (!G_random_element.inCoeffDomain())
1050 d0= totaldegree (G_random_element, Variable(2),
1051 Variable (G_random_element.level()));
1052 else
1053 d0= 0;
1054
1055 if (d0 == 0)
1056 {
1057 if (inextension)
1058 {
1059 ppA= GFMapDown (ppA, k);
1060 ppB= GFMapDown (ppB, k);
1061 setCharacteristic (p, k, gf_name_buf);
1062 }
1063 coF= N (ppA*(cA/gcdcAcB));
1064 coG= N (ppB*(cB/gcdcAcB));
1065 return N(gcdcAcB);
1066 }
1067 if (d0 > d)
1068 {
1069 if (!find (l, random_element))
1070 l.append (random_element);
1071 continue;
1072 }
1073
1074 G_random_element=
1075 (gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) *
1076 G_random_element;
1077
1078 coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
1079 *coF_random_element;
1080 coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
1081 *coG_random_element;
1082
1083 if (!G_random_element.inCoeffDomain())
1084 d0= totaldegree (G_random_element, Variable(2),
1085 Variable (G_random_element.level()));
1086 else
1087 d0= 0;
1088
1089 if (d0 < d)
1090 {
1091 m= gcdlcAlcB;
1092 newtonPoly= 1;
1093 G_m= 0;
1094 d= d0;
1095 coF_m= 0;
1096 coG_m= 0;
1097 }
1098
1099 TIMING_START (newton_interpolation);
1100 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1101 coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
1102 coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
1103 TIMING_END_AND_PRINT (newton_interpolation,
1104 "time for newton interpolation: ");
1105
1106 //termination test
1107 if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
1108 {
1109 TIMING_START (termination_test);
1110 if (gcdlcAlcB.isOne())
1111 cH= 1;
1112 else
1113 cH= uni_content (H);
1114 ppH= H/cH;
1115 ppH /= Lc (ppH);
1116 CanonicalForm lcppH= gcdlcAlcB/cH;
1117 CanonicalForm ccoF= lcppH/Lc (lcppH);
1118 CanonicalForm ccoG= lcppH/Lc (lcppH);
1119 ppCoF= coF/ccoF;
1120 ppCoG= coG/ccoG;
1121 if (inextension)
1122 {
1123 if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1124 (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1125 terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1126 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1127 {
1128 DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
1129 ppH= GFMapDown (ppH, k);
1130 ppCoF= GFMapDown (ppCoF, k);
1131 ppCoG= GFMapDown (ppCoG, k);
1132 DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
1133 coF= N ((cA/gcdcAcB)*ppCoF);
1134 coG= N ((cB/gcdcAcB)*ppCoG);
1135 setCharacteristic (p, k, gf_name_buf);
1136 TIMING_END_AND_PRINT (termination_test,
1137 "time for successful termination GF: ");
1138 return N(gcdcAcB*ppH);
1139 }
1140 }
1141 else
1142 {
1143 if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1144 (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1145 terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1146 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1147 {
1148 coF= N ((cA/gcdcAcB)*ppCoF);
1149 coG= N ((cB/gcdcAcB)*ppCoG);
1150 TIMING_END_AND_PRINT (termination_test,
1151 "time for successful termination GF: ");
1152 return N(gcdcAcB*ppH);
1153 }
1154 }
1155 TIMING_END_AND_PRINT (termination_test,
1156 "time for unsuccessful termination GF: ");
1157 }
1158
1159 G_m= H;
1160 coF_m= coF;
1161 coG_m= coG;
1162 newtonPoly= newtonPoly*(x - random_element);
1163 m= m*(x - random_element);
1164 if (!find (l, random_element))
1165 l.append (random_element);
1166 } while (1);
1167}
1168
1169static inline
1171FpRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
1172{
1173 fail= false;
1174 Variable x= F.mvar();
1175 FFRandom genFF;
1176 CanonicalForm random;
1177 int p= getCharacteristic();
1178 int bound= p;
1179 do
1180 {
1181 if (list.length() == bound)
1182 {
1183 fail= true;
1184 break;
1185 }
1186 if (list.length() < 1)
1187 random= 0;
1188 else
1189 {
1190 random= genFF.generate();
1191 while (find (list, random))
1192 random= genFF.generate();
1193 }
1194 if (F (random, x) == 0)
1195 {
1196 list.append (random);
1197 continue;
1198 }
1199 } while (find (list, random));
1200 return random;
1201}
1202
1203#if defined(HAVE_NTL) || defined(HAVE_FLINT)
1205modGCDFp (const CanonicalForm& F, const CanonicalForm& G,
1207 bool& topLevel, CFList& l);
1208#endif
1209
1210#if defined(HAVE_NTL) || defined(HAVE_FLINT)
1213 bool& topLevel, CFList& l)
1214{
1215 CanonicalForm dummy1, dummy2;
1216 CanonicalForm result= modGCDFp (F, G, dummy1, dummy2, topLevel, l);
1217 return result;
1218}
1219#endif
1220
1221#if defined(HAVE_NTL) || defined(HAVE_FLINT)
1225 bool& topLevel, CFList& l)
1226{
1227 CanonicalForm A= F;
1228 CanonicalForm B= G;
1229 if (F.isZero() && degree(G) > 0)
1230 {
1231 coF= 0;
1232 coG= Lc (G);
1233 return G/Lc(G);
1234 }
1235 else if (G.isZero() && degree (F) > 0)
1236 {
1237 coF= Lc (F);
1238 coG= 0;
1239 return F/Lc(F);
1240 }
1241 else if (F.isZero() && G.isZero())
1242 {
1243 coF= coG= 0;
1244 return F.genOne();
1245 }
1246 if (F.inBaseDomain() || G.inBaseDomain())
1247 {
1248 coF= F;
1249 coG= G;
1250 return F.genOne();
1251 }
1252 if (F.isUnivariate() && fdivides(F, G, coG))
1253 {
1254 coF= Lc (F);
1255 return F/Lc(F);
1256 }
1257 if (G.isUnivariate() && fdivides(G, F, coF))
1258 {
1259 coG= Lc (G);
1260 return G/Lc(G);
1261 }
1262 if (F == G)
1263 {
1264 coF= coG= Lc (F);
1265 return F/Lc(F);
1266 }
1267 CFMap M,N;
1268 int best_level= myCompress (A, B, M, N, topLevel);
1269
1270 if (best_level == 0)
1271 {
1272 coF= F;
1273 coG= G;
1274 return B.genOne();
1275 }
1276
1277 A= M(A);
1278 B= M(B);
1279
1280 Variable x= Variable (1);
1281
1282 //univariate case
1283 if (A.isUnivariate() && B.isUnivariate())
1284 {
1286 coF= N (A/result);
1287 coG= N (B/result);
1288 return N (result);
1289 }
1290
1291 CanonicalForm cA, cB; // content of A and B
1292 CanonicalForm ppA, ppB; // primitive part of A and B
1293 CanonicalForm gcdcAcB;
1294
1295 cA = uni_content (A);
1296 cB = uni_content (B);
1297 gcdcAcB= gcd (cA, cB);
1298 ppA= A/cA;
1299 ppB= B/cB;
1300
1301 CanonicalForm lcA, lcB; // leading coefficients of A and B
1302 CanonicalForm gcdlcAlcB;
1303 lcA= uni_lcoeff (ppA);
1304 lcB= uni_lcoeff (ppB);
1305
1306 gcdlcAlcB= gcd (lcA, lcB);
1307
1308 int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
1309 int d0;
1310
1311 if (d == 0)
1312 {
1313 coF= N (ppA*(cA/gcdcAcB));
1314 coG= N (ppB*(cB/gcdcAcB));
1315 return N(gcdcAcB);
1316 }
1317
1318 d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
1319
1320 if (d0 < d)
1321 d= d0;
1322
1323 if (d == 0)
1324 {
1325 coF= N (ppA*(cA/gcdcAcB));
1326 coG= N (ppB*(cB/gcdcAcB));
1327 return N(gcdcAcB);
1328 }
1329
1330 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
1331 CanonicalForm newtonPoly, coF_random_element, coG_random_element,
1332 coF_m, coG_m, ppCoF, ppCoG;
1333
1334 newtonPoly= 1;
1335 m= gcdlcAlcB;
1336 H= 0;
1337 coF= 0;
1338 coG= 0;
1339 G_m= 0;
1340 Variable alpha, V_buf, cleanUp;
1341 bool fail= false;
1342 bool inextension= false;
1343 topLevel= false;
1344 CFList source, dest;
1345 int bound1= degree (ppA, 1);
1346 int bound2= degree (ppB, 1);
1347 do
1348 {
1349 if (inextension)
1350 random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
1351 else
1352 random_element= FpRandomElement (m*lcA*lcB, l, fail);
1353
1354 if (!fail && !inextension)
1355 {
1356 CFList list;
1357 TIMING_START (gcd_recursion);
1358 G_random_element=
1359 modGCDFp (ppA (random_element,x), ppB (random_element,x),
1360 coF_random_element, coG_random_element, topLevel,
1361 list);
1362 TIMING_END_AND_PRINT (gcd_recursion,
1363 "time for recursive call: ");
1364 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1365 }
1366 else if (!fail && inextension)
1367 {
1368 CFList list;
1369 TIMING_START (gcd_recursion);
1370 G_random_element=
1371 modGCDFq (ppA (random_element, x), ppB (random_element, x),
1372 coF_random_element, coG_random_element, V_buf,
1373 list, topLevel);
1374 TIMING_END_AND_PRINT (gcd_recursion,
1375 "time for recursive call: ");
1376 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1377 }
1378 else if (fail && !inextension)
1379 {
1380 source= CFList();
1381 dest= CFList();
1382 CFList list;
1384 int deg= 2;
1385 bool initialized= false;
1386 do
1387 {
1388 mipo= randomIrredpoly (deg, x);
1389 if (initialized)
1390 setMipo (alpha, mipo);
1391 else
1392 alpha= rootOf (mipo);
1393 inextension= true;
1394 initialized= true;
1395 fail= false;
1396 random_element= randomElement (m*lcA*lcB, alpha, l, fail);
1397 deg++;
1398 } while (fail);
1399 list= CFList();
1400 V_buf= alpha;
1401 cleanUp= alpha;
1402 TIMING_START (gcd_recursion);
1403 G_random_element=
1404 modGCDFq (ppA (random_element, x), ppB (random_element, x),
1405 coF_random_element, coG_random_element, alpha,
1406 list, topLevel);
1407 TIMING_END_AND_PRINT (gcd_recursion,
1408 "time for recursive call: ");
1409 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1410 }
1411 else if (fail && inextension)
1412 {
1413 source= CFList();
1414 dest= CFList();
1415
1416 Variable V_buf3= V_buf;
1417 V_buf= chooseExtension (V_buf);
1418 bool prim_fail= false;
1419 Variable V_buf2;
1420 CanonicalForm prim_elem, im_prim_elem;
1421 prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
1422
1423 if (V_buf3 != alpha)
1424 {
1425 m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
1426 G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest);
1427 coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest);
1428 coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest);
1429 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
1430 source, dest);
1431 ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
1432 ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
1433 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
1434 dest);
1435 lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest);
1436 lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest);
1437 for (CFListIterator i= l; i.hasItem(); i++)
1438 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
1439 source, dest);
1440 }
1441
1442 ASSERT (!prim_fail, "failure in integer factorizer");
1443 if (prim_fail)
1444 ; //ERROR
1445 else
1446 im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
1447
1448 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
1449 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
1450
1451 for (CFListIterator i= l; i.hasItem(); i++)
1452 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
1453 im_prim_elem, source, dest);
1454 m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1455 G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1456 coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1457 coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1458 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
1459 source, dest);
1460 ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1461 ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1462 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
1463 source, dest);
1464 lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1465 lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1466 fail= false;
1467 random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
1468 DEBOUTLN (cerr, "fail= " << fail);
1469 CFList list;
1470 TIMING_START (gcd_recursion);
1471 G_random_element=
1472 modGCDFq (ppA (random_element, x), ppB (random_element, x),
1473 coF_random_element, coG_random_element, V_buf,
1474 list, topLevel);
1475 TIMING_END_AND_PRINT (gcd_recursion,
1476 "time for recursive call: ");
1477 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1478 alpha= V_buf;
1479 }
1480
1481 if (!G_random_element.inCoeffDomain())
1482 d0= totaldegree (G_random_element, Variable(2),
1483 Variable (G_random_element.level()));
1484 else
1485 d0= 0;
1486
1487 if (d0 == 0)
1488 {
1489 if (inextension)
1490 prune (cleanUp);
1491 coF= N (ppA*(cA/gcdcAcB));
1492 coG= N (ppB*(cB/gcdcAcB));
1493 return N(gcdcAcB);
1494 }
1495
1496 if (d0 > d)
1497 {
1498 if (!find (l, random_element))
1499 l.append (random_element);
1500 continue;
1501 }
1502
1503 G_random_element= (gcdlcAlcB(random_element,x)/uni_lcoeff(G_random_element))
1504 *G_random_element;
1505
1506 coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
1507 *coF_random_element;
1508 coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
1509 *coG_random_element;
1510
1511 if (!G_random_element.inCoeffDomain())
1512 d0= totaldegree (G_random_element, Variable(2),
1513 Variable (G_random_element.level()));
1514 else
1515 d0= 0;
1516
1517 if (d0 < d)
1518 {
1519 m= gcdlcAlcB;
1520 newtonPoly= 1;
1521 G_m= 0;
1522 d= d0;
1523 coF_m= 0;
1524 coG_m= 0;
1525 }
1526
1527 TIMING_START (newton_interpolation);
1528 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1529 coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
1530 coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
1531 TIMING_END_AND_PRINT (newton_interpolation,
1532 "time for newton_interpolation: ");
1533
1534 //termination test
1535 if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
1536 {
1537 TIMING_START (termination_test);
1538 if (gcdlcAlcB.isOne())
1539 cH= 1;
1540 else
1541 cH= uni_content (H);
1542 ppH= H/cH;
1543 ppH /= Lc (ppH);
1544 CanonicalForm lcppH= gcdlcAlcB/cH;
1545 CanonicalForm ccoF= lcppH/Lc (lcppH);
1546 CanonicalForm ccoG= lcppH/Lc (lcppH);
1547 ppCoF= coF/ccoF;
1548 ppCoG= coG/ccoG;
1549 DEBOUTLN (cerr, "ppH= " << ppH);
1550 if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1551 (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1552 terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1553 (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1554 {
1555 if (inextension)
1556 prune (cleanUp);
1557 coF= N ((cA/gcdcAcB)*ppCoF);
1558 coG= N ((cB/gcdcAcB)*ppCoG);
1559 TIMING_END_AND_PRINT (termination_test,
1560 "time for successful termination Fp: ");
1561 return N(gcdcAcB*ppH);
1562 }
1563 TIMING_END_AND_PRINT (termination_test,
1564 "time for unsuccessful termination Fp: ");
1565 }
1566
1567 G_m= H;
1568 coF_m= coF;
1569 coG_m= coG;
1570 newtonPoly= newtonPoly*(x - random_element);
1571 m= m*(x - random_element);
1572 if (!find (l, random_element))
1573 l.append (random_element);
1574 } while (1);
1575}
1576#endif
1577
1578CFArray
1580{
1581 int r= M.size();
1582 ASSERT (A.size() == r, "vector does not have right size");
1583
1584 if (r == 1)
1585 {
1586 CFArray result= CFArray (1);
1587 result [0]= A [0] / M [0];
1588 return result;
1589 }
1590 // check solvability
1591 bool notDistinct= false;
1592 for (int i= 0; i < r - 1; i++)
1593 {
1594 for (int j= i + 1; j < r; j++)
1595 {
1596 if (M [i] == M [j])
1597 {
1598 notDistinct= true;
1599 break;
1600 }
1601 }
1602 }
1603 if (notDistinct)
1604 return CFArray();
1605
1606 CanonicalForm master= 1;
1607 Variable x= Variable (1);
1608 for (int i= 0; i < r; i++)
1609 master *= x - M [i];
1610 CFList Pj;
1611 CanonicalForm tmp;
1612 for (int i= 0; i < r; i++)
1613 {
1614 tmp= master/(x - M [i]);
1615 tmp /= tmp (M [i], 1);
1616 Pj.append (tmp);
1617 }
1618 CFArray result= CFArray (r);
1619
1620 CFListIterator j= Pj;
1621 for (int i= 1; i <= r; i++, j++)
1622 {
1623 tmp= 0;
1624 for (int l= 0; l < A.size(); l++)
1625 tmp += A[l]*j.getItem()[l];
1626 result[i - 1]= tmp;
1627 }
1628 return result;
1629}
1630
1631CFArray
1633{
1634 int r= M.size();
1635 ASSERT (A.size() == r, "vector does not have right size");
1636 if (r == 1)
1637 {
1638 CFArray result= CFArray (1);
1639 result [0]= A[0] / M [0];
1640 return result;
1641 }
1642 // check solvability
1643 bool notDistinct= false;
1644 for (int i= 0; i < r - 1; i++)
1645 {
1646 for (int j= i + 1; j < r; j++)
1647 {
1648 if (M [i] == M [j])
1649 {
1650 notDistinct= true;
1651 break;
1652 }
1653 }
1654 }
1655 if (notDistinct)
1656 return CFArray();
1657
1658 CanonicalForm master= 1;
1659 Variable x= Variable (1);
1660 for (int i= 0; i < r; i++)
1661 master *= x - M [i];
1662 master *= x;
1663 CFList Pj;
1664 CanonicalForm tmp;
1665 for (int i= 0; i < r; i++)
1666 {
1667 tmp= master/(x - M [i]);
1668 tmp /= tmp (M [i], 1);
1669 Pj.append (tmp);
1670 }
1671
1672 CFArray result= CFArray (r);
1673
1674 CFListIterator j= Pj;
1675 for (int i= 1; i <= r; i++, j++)
1676 {
1677 tmp= 0;
1678
1679 for (int l= 1; l <= A.size(); l++)
1680 tmp += A[l - 1]*j.getItem()[l];
1681 result[i - 1]= tmp;
1682 }
1683 return result;
1684}
1685
1686/// M in row echolon form, rk rank of M
1687CFArray
1688readOffSolution (const CFMatrix& M, const long rk)
1689{
1690 CFArray result= CFArray (rk);
1691 CanonicalForm tmp1, tmp2, tmp3;
1692 for (int i= rk; i >= 1; i--)
1693 {
1694 tmp3= 0;
1695 tmp1= M (i, M.columns());
1696 for (int j= M.columns() - 1; j >= 1; j--)
1697 {
1698 tmp2= M (i, j);
1699 if (j == i)
1700 break;
1701 else
1702 tmp3 += tmp2*result[j - 1];
1703 }
1704 result[i - 1]= (tmp1 - tmp3)/tmp2;
1705 }
1706 return result;
1707}
1708
1709CFArray
1710readOffSolution (const CFMatrix& M, const CFArray& L, const CFArray& partialSol)
1711{
1712 CFArray result= CFArray (M.rows());
1713 CanonicalForm tmp1, tmp2, tmp3;
1714 int k;
1715 for (int i= M.rows(); i >= 1; i--)
1716 {
1717 tmp3= 0;
1718 tmp1= L[i - 1];
1719 k= 0;
1720 for (int j= M.columns(); j >= 1; j--, k++)
1721 {
1722 tmp2= M (i, j);
1723 if (j == i)
1724 break;
1725 else
1726 {
1727 if (k > partialSol.size() - 1)
1728 tmp3 += tmp2*result[j - 1];
1729 else
1730 tmp3 += tmp2*partialSol[partialSol.size() - k - 1];
1731 }
1732 }
1733 result [i - 1]= (tmp1 - tmp3)/tmp2;
1734 }
1735 return result;
1736}
1737
1738long
1740{
1741 ASSERT (L.size() <= M.rows(), "dimension exceeded");
1742 CFMatrix *N;
1743 N= new CFMatrix (M.rows(), M.columns() + 1);
1744
1745 for (int i= 1; i <= M.rows(); i++)
1746 for (int j= 1; j <= M.columns(); j++)
1747 (*N) (i, j)= M (i, j);
1748
1749 int j= 1;
1750 for (int i= 0; i < L.size(); i++, j++)
1751 (*N) (j, M.columns() + 1)= L[i];
1752#ifdef HAVE_FLINT
1753 nmod_mat_t FLINTN;
1755 long rk= nmod_mat_rref (FLINTN);
1756
1757 delete N;
1759 nmod_mat_clear (FLINTN);
1760#else
1761 int p= getCharacteristic ();
1762 if (fac_NTL_char != p)
1763 {
1764 fac_NTL_char= p;
1765 zz_p::init (p);
1766 }
1767 mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1768 delete N;
1769 long rk= gauss (*NTLN);
1770
1772 delete NTLN;
1773#endif
1774
1775 L= CFArray (M.rows());
1776 for (int i= 0; i < M.rows(); i++)
1777 L[i]= (*N) (i + 1, M.columns() + 1);
1778 M= (*N) (1, M.rows(), 1, M.columns());
1779 delete N;
1780 return rk;
1781}
1782
1783long
1785{
1786 ASSERT (L.size() <= M.rows(), "dimension exceeded");
1787 CFMatrix *N;
1788 N= new CFMatrix (M.rows(), M.columns() + 1);
1789
1790 for (int i= 1; i <= M.rows(); i++)
1791 for (int j= 1; j <= M.columns(); j++)
1792 (*N) (i, j)= M (i, j);
1793
1794 int j= 1;
1795 for (int i= 0; i < L.size(); i++, j++)
1796 (*N) (j, M.columns() + 1)= L[i];
1797 #ifdef HAVE_FLINT
1798 // convert mipo
1799 nmod_poly_t mipo1;
1801 fq_nmod_ctx_t ctx;
1802 fq_nmod_ctx_init_modulus(ctx,mipo1,"t");
1803 nmod_poly_clear(mipo1);
1804 // convert matrix
1805 fq_nmod_mat_t FLINTN;
1806 convertFacCFMatrix2Fq_nmod_mat_t (FLINTN, ctx, *N);
1807 // rank
1808 long rk= fq_nmod_mat_rref (FLINTN,ctx);
1809 // clean up
1810 fq_nmod_mat_clear (FLINTN,ctx);
1811 fq_nmod_ctx_clear(ctx);
1812 #elif defined(HAVE_NTL)
1813 int p= getCharacteristic ();
1814 if (fac_NTL_char != p)
1815 {
1816 fac_NTL_char= p;
1817 zz_p::init (p);
1818 }
1819 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1820 zz_pE::init (NTLMipo);
1821 mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1822 long rk= gauss (*NTLN);
1824 delete NTLN;
1825 #else
1826 factoryError("NTL/FLINT missing: gaussianElimFq");
1827 #endif
1828
1829 M= (*N) (1, M.rows(), 1, M.columns());
1830 L= CFArray (M.rows());
1831 for (int i= 0; i < M.rows(); i++)
1832 L[i]= (*N) (i + 1, M.columns() + 1);
1833
1834 delete N;
1835 return rk;
1836}
1837
1838CFArray
1840{
1841 ASSERT (L.size() <= M.rows(), "dimension exceeded");
1842 CFMatrix *N;
1843 N= new CFMatrix (M.rows(), M.columns() + 1);
1844
1845 for (int i= 1; i <= M.rows(); i++)
1846 for (int j= 1; j <= M.columns(); j++)
1847 (*N) (i, j)= M (i, j);
1848
1849 int j= 1;
1850 for (int i= 0; i < L.size(); i++, j++)
1851 (*N) (j, M.columns() + 1)= L[i];
1852
1853#ifdef HAVE_FLINT
1854 nmod_mat_t FLINTN;
1856 long rk= nmod_mat_rref (FLINTN);
1857#else
1858 int p= getCharacteristic ();
1859 if (fac_NTL_char != p)
1860 {
1861 fac_NTL_char= p;
1862 zz_p::init (p);
1863 }
1864 mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1865 long rk= gauss (*NTLN);
1866#endif
1867 delete N;
1868 if (rk != M.columns())
1869 {
1870#ifdef HAVE_FLINT
1871 nmod_mat_clear (FLINTN);
1872#else
1873 delete NTLN;
1874#endif
1875 return CFArray();
1876 }
1877#ifdef HAVE_FLINT
1879 nmod_mat_clear (FLINTN);
1880#else
1882 delete NTLN;
1883#endif
1884 CFArray A= readOffSolution (*N, rk);
1885
1886 delete N;
1887 return A;
1888}
1889
1890CFArray
1891solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)
1892{
1893 ASSERT (L.size() <= M.rows(), "dimension exceeded");
1894 CFMatrix *N;
1895 N= new CFMatrix (M.rows(), M.columns() + 1);
1896
1897 for (int i= 1; i <= M.rows(); i++)
1898 for (int j= 1; j <= M.columns(); j++)
1899 (*N) (i, j)= M (i, j);
1900 int j= 1;
1901 for (int i= 0; i < L.size(); i++, j++)
1902 (*N) (j, M.columns() + 1)= L[i];
1903 #ifdef HAVE_FLINT
1904 // convert mipo
1905 nmod_poly_t mipo1;
1907 fq_nmod_ctx_t ctx;
1908 fq_nmod_ctx_init_modulus(ctx,mipo1,"t");
1909 nmod_poly_clear(mipo1);
1910 // convert matrix
1911 fq_nmod_mat_t FLINTN;
1912 convertFacCFMatrix2Fq_nmod_mat_t (FLINTN, ctx, *N);
1913 // rank
1914 long rk= fq_nmod_mat_rref (FLINTN,ctx);
1915 #elif defined(HAVE_NTL)
1916 int p= getCharacteristic ();
1917 if (fac_NTL_char != p)
1918 {
1919 fac_NTL_char= p;
1920 zz_p::init (p);
1921 }
1922 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1923 zz_pE::init (NTLMipo);
1924 mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1925 long rk= gauss (*NTLN);
1926 #else
1927 factoryError("NTL/FLINT missing: solveSystemFq");
1928 #endif
1929
1930 delete N;
1931 if (rk != M.columns())
1932 {
1933 #if defined(HAVE_NTL) && !defined(HAVE_FLINT)
1934 delete NTLN;
1935 #endif
1936 return CFArray();
1937 }
1938 #ifdef HAVE_FLINT
1939 // convert and clean up
1941 fq_nmod_mat_clear (FLINTN,ctx);
1942 fq_nmod_ctx_clear(ctx);
1943 #elif defined(HAVE_NTL)
1945 delete NTLN;
1946 #endif
1947
1948 CFArray A= readOffSolution (*N, rk);
1949
1950 delete N;
1951 return A;
1952}
1953#endif
1954
1955CFArray
1957{
1958 if (F.inCoeffDomain())
1959 {
1960 CFArray result= CFArray (1);
1961 result [0]= 1;
1962 return result;
1963 }
1964 if (F.isUnivariate())
1965 {
1966 CFArray result= CFArray (size(F));
1967 int j= 0;
1968 for (CFIterator i= F; i.hasTerms(); i++, j++)
1969 result[j]= power (F.mvar(), i.exp());
1970 return result;
1971 }
1972 int numMon= size (F);
1973 CFArray result= CFArray (numMon);
1974 int j= 0;
1975 CFArray recResult;
1976 Variable x= F.mvar();
1977 CanonicalForm powX;
1978 for (CFIterator i= F; i.hasTerms(); i++)
1979 {
1980 powX= power (x, i.exp());
1981 recResult= getMonoms (i.coeff());
1982 for (int k= 0; k < recResult.size(); k++)
1983 result[j+k]= powX*recResult[k];
1984 j += recResult.size();
1985 }
1986 return result;
1987}
1988
1989#if defined(HAVE_NTL) || defined(HAVE_FLINT)
1990CFArray
1992{
1993 if (F.inCoeffDomain())
1994 {
1995 CFArray result= CFArray (1);
1996 result [0]= F;
1997 return result;
1998 }
1999 if (F.isUnivariate())
2000 {
2001 ASSERT (evalPoints.length() == 1,
2002 "expected an eval point with only one component");
2003 CFArray result= CFArray (size(F));
2004 int j= 0;
2006 for (CFIterator i= F; i.hasTerms(); i++, j++)
2007 result[j]= power (evalPoint, i.exp());
2008 return result;
2009 }
2010 int numMon= size (F);
2011 CFArray result= CFArray (numMon);
2012 int j= 0;
2015 buf.removeLast();
2016 CFArray recResult;
2017 CanonicalForm powEvalPoint;
2018 for (CFIterator i= F; i.hasTerms(); i++)
2019 {
2020 powEvalPoint= power (evalPoint, i.exp());
2021 recResult= evaluateMonom (i.coeff(), buf);
2022 for (int k= 0; k < recResult.size(); k++)
2023 result[j+k]= powEvalPoint*recResult[k];
2024 j += recResult.size();
2025 }
2026 return result;
2027}
2028
2029CFArray
2031{
2032 CFArray result= A.size();
2033 CanonicalForm tmp;
2034 int k;
2035 for (int i= 0; i < A.size(); i++)
2036 {
2037 tmp= A[i];
2038 k= 1;
2039 for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
2040 tmp= tmp (j.getItem(), k);
2041 result[i]= tmp;
2042 }
2043 return result;
2044}
2045
2046CFList
2049 const CanonicalForm& LCF, const bool& GF,
2050 const Variable& alpha, bool& fail, CFList& list
2051 )
2052{
2053 int k= tmax (F.level(), G.level()) - 1;
2054 Variable x= Variable (1);
2055 CFList result;
2056 FFRandom genFF;
2057 GFRandom genGF;
2058 int p= getCharacteristic ();
2059 double bound;
2060 if (alpha != Variable (1))
2061 {
2062 bound= pow ((double) p, (double) degree (getMipo(alpha)));
2063 bound= pow (bound, (double) k);
2064 }
2065 else if (GF)
2066 {
2067 bound= pow ((double) p, (double) getGFDegree());
2068 bound= pow ((double) bound, (double) k);
2069 }
2070 else
2071 bound= pow ((double) p, (double) k);
2072
2073 CanonicalForm random;
2074 int j;
2075 bool zeroOneOccured= false;
2076 bool allEqual= false;
2078 do
2079 {
2080 random= 0;
2081 // possible overflow if list.length() does not fit into a int
2082 if (list.length() >= bound)
2083 {
2084 fail= true;
2085 break;
2086 }
2087 for (int i= 0; i < k; i++)
2088 {
2089 if (GF)
2090 {
2091 result.append (genGF.generate());
2092 random += result.getLast()*power (x, i);
2093 }
2094 else if (alpha.level() != 1)
2095 {
2096 AlgExtRandomF genAlgExt (alpha);
2097 result.append (genAlgExt.generate());
2098 random += result.getLast()*power (x, i);
2099 }
2100 else
2101 {
2102 result.append (genFF.generate());
2103 random += result.getLast()*power (x, i);
2104 }
2105 if (result.getLast().isOne() || result.getLast().isZero())
2106 zeroOneOccured= true;
2107 }
2108 if (find (list, random))
2109 {
2110 zeroOneOccured= false;
2111 allEqual= false;
2112 result= CFList();
2113 continue;
2114 }
2115 if (zeroOneOccured)
2116 {
2117 list.append (random);
2118 zeroOneOccured= false;
2119 allEqual= false;
2120 result= CFList();
2121 continue;
2122 }
2123 // no zero at this point
2124 if (k > 1)
2125 {
2126 allEqual= true;
2127 CFIterator iter= random;
2128 buf= iter.coeff();
2129 iter++;
2130 for (; iter.hasTerms(); iter++)
2131 if (buf != iter.coeff())
2132 allEqual= false;
2133 }
2134 if (allEqual)
2135 {
2136 list.append (random);
2137 allEqual= false;
2138 zeroOneOccured= false;
2139 result= CFList();
2140 continue;
2141 }
2142
2143 Feval= F;
2144 Geval= G;
2145 CanonicalForm LCeval= LCF;
2146 j= 1;
2147 for (CFListIterator i= result; i.hasItem(); i++, j++)
2148 {
2149 Feval= Feval (i.getItem(), j);
2150 Geval= Geval (i.getItem(), j);
2151 LCeval= LCeval (i.getItem(), j);
2152 }
2153
2154 if (LCeval.isZero())
2155 {
2156 if (!find (list, random))
2157 list.append (random);
2158 zeroOneOccured= false;
2159 allEqual= false;
2160 result= CFList();
2161 continue;
2162 }
2163
2164 if (list.length() >= bound)
2165 {
2166 fail= true;
2167 break;
2168 }
2169 } while (find (list, random));
2170
2171 return result;
2172}
2173
2174/// multiply two lists componentwise
2175void mult (CFList& L1, const CFList& L2)
2176{
2177 ASSERT (L1.length() == L2.length(), "lists of the same size expected");
2178
2179 CFListIterator j= L2;
2180 for (CFListIterator i= L1; i.hasItem(); i++, j++)
2181 i.getItem() *= j.getItem();
2182}
2183
2185 CanonicalForm& Beval, const CFList& L)
2186{
2187 Aeval= A;
2188 Beval= B;
2189 int j= 1;
2190 for (CFListIterator i= L; i.hasItem(); i++, j++)
2191 {
2192 Aeval= Aeval (i.getItem(), j);
2193 Beval= Beval (i.getItem(), j);
2194 }
2195}
2196
2199 const CanonicalForm& skeleton, const Variable& alpha,
2200 bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2201 )
2202{
2203 CanonicalForm A= F;
2204 CanonicalForm B= G;
2205 if (F.isZero() && degree(G) > 0) return G/Lc(G);
2206 else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2207 else if (F.isZero() && G.isZero()) return F.genOne();
2208 if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2209 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2210 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2211 if (F == G) return F/Lc(F);
2212
2213 ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2214 ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2215
2216 CFMap M,N;
2217 int best_level= myCompress (A, B, M, N, false);
2218
2219 if (best_level == 0)
2220 return B.genOne();
2221
2222 A= M(A);
2223 B= M(B);
2224
2225 Variable x= Variable (1);
2226
2227 //univariate case
2228 if (A.isUnivariate() && B.isUnivariate())
2229 return N (gcd (A, B));
2230
2231 CanonicalForm skel= M(skeleton);
2232 CanonicalForm cA, cB; // content of A and B
2233 CanonicalForm ppA, ppB; // primitive part of A and B
2234 CanonicalForm gcdcAcB;
2235 cA = uni_content (A);
2236 cB = uni_content (B);
2237 gcdcAcB= gcd (cA, cB);
2238 ppA= A/cA;
2239 ppB= B/cB;
2240
2241 CanonicalForm lcA, lcB; // leading coefficients of A and B
2242 CanonicalForm gcdlcAlcB;
2243 lcA= uni_lcoeff (ppA);
2244 lcB= uni_lcoeff (ppB);
2245
2246 if (fdivides (lcA, lcB))
2247 {
2248 if (fdivides (A, B))
2249 return F/Lc(F);
2250 }
2251 if (fdivides (lcB, lcA))
2252 {
2253 if (fdivides (B, A))
2254 return G/Lc(G);
2255 }
2256
2257 gcdlcAlcB= gcd (lcA, lcB);
2258 int skelSize= size (skel, skel.mvar());
2259
2260 int j= 0;
2261 int biggestSize= 0;
2262
2263 for (CFIterator i= skel; i.hasTerms(); i++, j++)
2264 biggestSize= tmax (biggestSize, size (i.coeff()));
2265
2266 CanonicalForm g, Aeval, Beval;
2267
2269 bool evalFail= false;
2270 CFList list;
2271 bool GF= false;
2272 CanonicalForm LCA= LC (A);
2273 CanonicalForm tmp;
2274 CFArray gcds= CFArray (biggestSize);
2275 CFList * pEvalPoints= new CFList [biggestSize];
2276 Variable V_buf= alpha, V_buf4= alpha;
2277 CFList source, dest;
2278 CanonicalForm prim_elem, im_prim_elem;
2279 CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2280 for (int i= 0; i < biggestSize; i++)
2281 {
2282 if (i == 0)
2283 evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
2284 list);
2285 else
2286 {
2287 mult (evalPoints, pEvalPoints [0]);
2288 eval (A, B, Aeval, Beval, evalPoints);
2289 }
2290
2291 if (evalFail)
2292 {
2293 if (V_buf.level() != 1)
2294 {
2295 do
2296 {
2297 Variable V_buf3= V_buf;
2298 V_buf= chooseExtension (V_buf);
2299 source= CFList();
2300 dest= CFList();
2301
2302 bool prim_fail= false;
2303 Variable V_buf2;
2304 prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2305 if (V_buf4 == alpha && alpha.level() != 1)
2306 prim_elem_alpha= prim_elem;
2307
2308 ASSERT (!prim_fail, "failure in integer factorizer");
2309 if (prim_fail)
2310 ; //ERROR
2311 else
2312 im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2313
2314 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2315 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2316
2317 if (V_buf4 == alpha && alpha.level() != 1)
2318 im_prim_elem_alpha= im_prim_elem;
2319 else if (alpha.level() != 1)
2320 im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2321 prim_elem, im_prim_elem, source, dest);
2322
2323 for (CFListIterator j= list; j.hasItem(); j++)
2324 j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2325 im_prim_elem, source, dest);
2326 for (int k= 0; k < i; k++)
2327 {
2328 for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
2329 j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2330 im_prim_elem, source, dest);
2331 gcds[k]= mapUp (gcds[k], V_buf4, V_buf, prim_elem, im_prim_elem,
2332 source, dest);
2333 }
2334
2335 if (alpha.level() != 1)
2336 {
2337 A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2338 B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2339 }
2340 V_buf4= V_buf;
2341 evalFail= false;
2342 evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2343 evalFail, list);
2344 } while (evalFail);
2345 }
2346 else
2347 {
2349 int deg= 2;
2350 bool initialized= false;
2351 do
2352 {
2353 mipo= randomIrredpoly (deg, x);
2354 if (initialized)
2355 setMipo (V_buf, mipo);
2356 else
2357 V_buf= rootOf (mipo);
2358 evalFail= false;
2359 initialized= true;
2360 evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2361 evalFail, list);
2362 deg++;
2363 } while (evalFail);
2364 V_buf4= V_buf;
2365 }
2366 }
2367
2368 g= gcd (Aeval, Beval);
2369 g /= Lc (g);
2370
2371 if (degree (g) != degree (skel) || (skelSize != size (g)))
2372 {
2373 delete[] pEvalPoints;
2374 fail= true;
2375 if (alpha.level() != 1 && V_buf != alpha)
2376 prune1 (alpha);
2377 return 0;
2378 }
2379 CFIterator l= skel;
2380 for (CFIterator k= g; k.hasTerms(); k++, l++)
2381 {
2382 if (k.exp() != l.exp())
2383 {
2384 delete[] pEvalPoints;
2385 fail= true;
2386 if (alpha.level() != 1 && V_buf != alpha)
2387 prune1 (alpha);
2388 return 0;
2389 }
2390 }
2391 pEvalPoints[i]= evalPoints;
2392 gcds[i]= g;
2393
2394 tmp= 0;
2395 int j= 0;
2396 for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2397 tmp += k.getItem()*power (x, j);
2398 list.append (tmp);
2399 }
2400
2401 if (Monoms.size() == 0)
2402 Monoms= getMonoms (skel);
2403
2404 coeffMonoms= new CFArray [skelSize];
2405 j= 0;
2406 for (CFIterator i= skel; i.hasTerms(); i++, j++)
2407 coeffMonoms[j]= getMonoms (i.coeff());
2408
2409 CFArray* pL= new CFArray [skelSize];
2410 CFArray* pM= new CFArray [skelSize];
2411 for (int i= 0; i < biggestSize; i++)
2412 {
2413 CFIterator l= gcds [i];
2414 evalPoints= pEvalPoints [i];
2415 for (int k= 0; k < skelSize; k++, l++)
2416 {
2417 if (i == 0)
2418 pL[k]= CFArray (biggestSize);
2419 pL[k] [i]= l.coeff();
2420
2421 if (i == 0)
2422 pM[k]= evaluate (coeffMonoms [k], evalPoints);
2423 }
2424 }
2425
2426 CFArray solution;
2428 int ind= 0;
2429 CFArray bufArray;
2430 CFMatrix Mat;
2431 for (int k= 0; k < skelSize; k++)
2432 {
2433 if (biggestSize != coeffMonoms[k].size())
2434 {
2435 bufArray= CFArray (coeffMonoms[k].size());
2436 for (int i= 0; i < coeffMonoms[k].size(); i++)
2437 bufArray [i]= pL[k] [i];
2438 solution= solveGeneralVandermonde (pM [k], bufArray);
2439 }
2440 else
2441 solution= solveGeneralVandermonde (pM [k], pL[k]);
2442
2443 if (solution.size() == 0)
2444 {
2445 delete[] pEvalPoints;
2446 delete[] pM;
2447 delete[] pL;
2448 delete[] coeffMonoms;
2449 fail= true;
2450 if (alpha.level() != 1 && V_buf != alpha)
2451 prune1 (alpha);
2452 return 0;
2453 }
2454 for (int l= 0; l < solution.size(); l++)
2455 result += solution[l]*Monoms [ind + l];
2456 ind += solution.size();
2457 }
2458
2459 delete[] pEvalPoints;
2460 delete[] pM;
2461 delete[] pL;
2462 delete[] coeffMonoms;
2463
2464 if (alpha.level() != 1 && V_buf != alpha)
2465 {
2466 CFList u, v;
2467 result= mapDown (result, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
2468 prune1 (alpha);
2469 }
2470
2471 result= N(result);
2472 if (fdivides (result, F) && fdivides (result, G))
2473 return result;
2474 else
2475 {
2476 fail= true;
2477 return 0;
2478 }
2479}
2480
2483 const CanonicalForm& skeleton, const Variable& alpha,
2484 bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2485 )
2486{
2487 CanonicalForm A= F;
2488 CanonicalForm B= G;
2489 if (F.isZero() && degree(G) > 0) return G/Lc(G);
2490 else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2491 else if (F.isZero() && G.isZero()) return F.genOne();
2492 if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2493 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2494 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2495 if (F == G) return F/Lc(F);
2496
2497 ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2498 ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2499
2500 CFMap M,N;
2501 int best_level= myCompress (A, B, M, N, false);
2502
2503 if (best_level == 0)
2504 return B.genOne();
2505
2506 A= M(A);
2507 B= M(B);
2508
2509 Variable x= Variable (1);
2510
2511 //univariate case
2512 if (A.isUnivariate() && B.isUnivariate())
2513 return N (gcd (A, B));
2514
2515 CanonicalForm skel= M(skeleton);
2516
2517 CanonicalForm cA, cB; // content of A and B
2518 CanonicalForm ppA, ppB; // primitive part of A and B
2519 CanonicalForm gcdcAcB;
2520 cA = uni_content (A);
2521 cB = uni_content (B);
2522 gcdcAcB= gcd (cA, cB);
2523 ppA= A/cA;
2524 ppB= B/cB;
2525
2526 CanonicalForm lcA, lcB; // leading coefficients of A and B
2527 CanonicalForm gcdlcAlcB;
2528 lcA= uni_lcoeff (ppA);
2529 lcB= uni_lcoeff (ppB);
2530
2531 if (fdivides (lcA, lcB))
2532 {
2533 if (fdivides (A, B))
2534 return F/Lc(F);
2535 }
2536 if (fdivides (lcB, lcA))
2537 {
2538 if (fdivides (B, A))
2539 return G/Lc(G);
2540 }
2541
2542 gcdlcAlcB= gcd (lcA, lcB);
2543 int skelSize= size (skel, skel.mvar());
2544
2545 int j= 0;
2546 int biggestSize= 0;
2547 int bufSize;
2548 int numberUni= 0;
2549 for (CFIterator i= skel; i.hasTerms(); i++, j++)
2550 {
2551 bufSize= size (i.coeff());
2552 biggestSize= tmax (biggestSize, bufSize);
2553 numberUni += bufSize;
2554 }
2555 numberUni--;
2556 numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
2557 biggestSize= tmax (biggestSize , numberUni);
2558
2559 numberUni= biggestSize + size (LC(skel)) - 1;
2560 int biggestSize2= tmax (numberUni, biggestSize);
2561
2562 CanonicalForm g, Aeval, Beval;
2563
2565 CFArray coeffEval;
2566 bool evalFail= false;
2567 CFList list;
2568 bool GF= false;
2569 CanonicalForm LCA= LC (A);
2570 CanonicalForm tmp;
2571 CFArray gcds= CFArray (biggestSize);
2572 CFList * pEvalPoints= new CFList [biggestSize];
2573 Variable V_buf= alpha, V_buf4= alpha;
2574 CFList source, dest;
2575 CanonicalForm prim_elem, im_prim_elem;
2576 CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2577 for (int i= 0; i < biggestSize; i++)
2578 {
2579 if (i == 0)
2580 {
2581 if (getCharacteristic() > 3)
2582 evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2583 evalFail, list);
2584 else
2585 evalFail= true;
2586
2587 if (evalFail)
2588 {
2589 if (V_buf.level() != 1)
2590 {
2591 do
2592 {
2593 Variable V_buf3= V_buf;
2594 V_buf= chooseExtension (V_buf);
2595 source= CFList();
2596 dest= CFList();
2597
2598 bool prim_fail= false;
2599 Variable V_buf2;
2600 prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2601 if (V_buf4 == alpha && alpha.level() != 1)
2602 prim_elem_alpha= prim_elem;
2603
2604 ASSERT (!prim_fail, "failure in integer factorizer");
2605 if (prim_fail)
2606 ; //ERROR
2607 else
2608 im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2609
2610 DEBOUTLN (cerr, "getMipo (V_buf)= " << getMipo (V_buf));
2611 DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
2612
2613 if (V_buf4 == alpha && alpha.level() != 1)
2614 im_prim_elem_alpha= im_prim_elem;
2615 else if (alpha.level() != 1)
2616 im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2617 prim_elem, im_prim_elem, source, dest);
2618
2619 for (CFListIterator i= list; i.hasItem(); i++)
2620 i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
2621 im_prim_elem, source, dest);
2622 if (alpha.level() != 1)
2623 {
2624 A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2625 B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2626 }
2627 evalFail= false;
2628 V_buf4= V_buf;
2629 evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2630 evalFail, list);
2631 } while (evalFail);
2632 }
2633 else
2634 {
2636 int deg= 2;
2637 bool initialized= false;
2638 do
2639 {
2640 mipo= randomIrredpoly (deg, x);
2641 if (initialized)
2642 setMipo (V_buf, mipo);
2643 else
2644 V_buf= rootOf (mipo);
2645 evalFail= false;
2646 initialized= true;
2647 evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2648 evalFail, list);
2649 deg++;
2650 } while (evalFail);
2651 V_buf4= V_buf;
2652 }
2653 }
2654 }
2655 else
2656 {
2657 mult (evalPoints, pEvalPoints[0]);
2658 eval (A, B, Aeval, Beval, evalPoints);
2659 }
2660
2661 g= gcd (Aeval, Beval);
2662 g /= Lc (g);
2663
2664 if (degree (g) != degree (skel) || (skelSize != size (g)))
2665 {
2666 delete[] pEvalPoints;
2667 fail= true;
2668 if (alpha.level() != 1 && V_buf != alpha)
2669 prune1 (alpha);
2670 return 0;
2671 }
2672 CFIterator l= skel;
2673 for (CFIterator k= g; k.hasTerms(); k++, l++)
2674 {
2675 if (k.exp() != l.exp())
2676 {
2677 delete[] pEvalPoints;
2678 fail= true;
2679 if (alpha.level() != 1 && V_buf != alpha)
2680 prune1 (alpha);
2681 return 0;
2682 }
2683 }
2684 pEvalPoints[i]= evalPoints;
2685 gcds[i]= g;
2686
2687 tmp= 0;
2688 int j= 0;
2689 for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2690 tmp += k.getItem()*power (x, j);
2691 list.append (tmp);
2692 }
2693
2694 if (Monoms.size() == 0)
2695 Monoms= getMonoms (skel);
2696
2697 coeffMonoms= new CFArray [skelSize];
2698
2699 j= 0;
2700 for (CFIterator i= skel; i.hasTerms(); i++, j++)
2701 coeffMonoms[j]= getMonoms (i.coeff());
2702
2703 int minimalColumnsIndex;
2704 if (skelSize > 1)
2705 minimalColumnsIndex= 1;
2706 else
2707 minimalColumnsIndex= 0;
2708 int minimalColumns=-1;
2709
2710 CFArray* pM= new CFArray [skelSize];
2711 CFMatrix Mat;
2712 // find the Matrix with minimal number of columns
2713 for (int i= 0; i < skelSize; i++)
2714 {
2715 pM[i]= CFArray (coeffMonoms[i].size());
2716 if (i == 1)
2717 minimalColumns= coeffMonoms[i].size();
2718 if (i > 1)
2719 {
2720 if (minimalColumns > coeffMonoms[i].size())
2721 {
2722 minimalColumns= coeffMonoms[i].size();
2723 minimalColumnsIndex= i;
2724 }
2725 }
2726 }
2727 CFMatrix* pMat= new CFMatrix [2];
2728 pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());
2729 pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());
2730 CFArray* pL= new CFArray [skelSize];
2731 for (int i= 0; i < biggestSize; i++)
2732 {
2733 CFIterator l= gcds [i];
2734 evalPoints= pEvalPoints [i];
2735 for (int k= 0; k < skelSize; k++, l++)
2736 {
2737 if (i == 0)
2738 pL[k]= CFArray (biggestSize);
2739 pL[k] [i]= l.coeff();
2740
2741 if (i == 0 && k != 0 && k != minimalColumnsIndex)
2742 pM[k]= evaluate (coeffMonoms[k], evalPoints);
2743 else if (i == 0 && (k == 0 || k == minimalColumnsIndex))
2744 coeffEval= evaluate (coeffMonoms[k], evalPoints);
2745 if (i > 0 && (k == 0 || k == minimalColumnsIndex))
2746 coeffEval= evaluate (coeffMonoms[k], evalPoints);
2747
2748 if (k == 0)
2749 {
2750 if (pMat[k].rows() >= i + 1)
2751 {
2752 for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2753 pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2754 }
2755 }
2756 if (k == minimalColumnsIndex)
2757 {
2758 if (pMat[1].rows() >= i + 1)
2759 {
2760 for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2761 pMat[1] (i + 1, ind)= coeffEval[ind - 1];
2762 }
2763 }
2764 }
2765 }
2766
2767 CFArray solution;
2769 int ind= 1;
2770 int matRows, matColumns;
2771 matRows= pMat[1].rows();
2772 matColumns= pMat[0].columns() - 1;
2773 matColumns += pMat[1].columns();
2774
2775 Mat= CFMatrix (matRows, matColumns);
2776 for (int i= 1; i <= matRows; i++)
2777 for (int j= 1; j <= pMat[1].columns(); j++)
2778 Mat (i, j)= pMat[1] (i, j);
2779
2780 for (int j= pMat[1].columns() + 1; j <= matColumns;
2781 j++, ind++)
2782 {
2783 for (int i= 1; i <= matRows; i++)
2784 Mat (i, j)= (-pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i - 1];
2785 }
2786
2787 CFArray firstColumn= CFArray (pMat[0].rows());
2788 for (int i= 0; i < pMat[0].rows(); i++)
2789 firstColumn [i]= pMat[0] (i + 1, 1);
2790
2791 CFArray bufArray= pL[minimalColumnsIndex];
2792
2793 for (int i= 0; i < biggestSize; i++)
2794 pL[minimalColumnsIndex] [i] *= firstColumn[i];
2795
2796 CFMatrix bufMat= pMat[1];
2797 pMat[1]= Mat;
2798
2799 if (V_buf.level() != 1)
2800 solution= solveSystemFq (pMat[1],
2801 pL[minimalColumnsIndex], V_buf);
2802 else
2803 solution= solveSystemFp (pMat[1],
2804 pL[minimalColumnsIndex]);
2805
2806 if (solution.size() == 0)
2807 { //Javadi's method failed, try deKleine, Monagan, Wittkopf instead
2808 CFMatrix bufMat0= pMat[0];
2809 delete [] pMat;
2810 pMat= new CFMatrix [skelSize];
2811 pL[minimalColumnsIndex]= bufArray;
2812 CFList* bufpEvalPoints= NULL;
2813 CFArray bufGcds;
2814 if (biggestSize != biggestSize2)
2815 {
2816 bufpEvalPoints= pEvalPoints;
2817 pEvalPoints= new CFList [biggestSize2];
2818 bufGcds= gcds;
2819 gcds= CFArray (biggestSize2);
2820 for (int i= 0; i < biggestSize; i++)
2821 {
2822 pEvalPoints[i]= bufpEvalPoints [i];
2823 gcds[i]= bufGcds[i];
2824 }
2825 for (int i= 0; i < biggestSize2 - biggestSize; i++)
2826 {
2827 mult (evalPoints, pEvalPoints[0]);
2828 eval (A, B, Aeval, Beval, evalPoints);
2829 g= gcd (Aeval, Beval);
2830 g /= Lc (g);
2831 if (degree (g) != degree (skel) || (skelSize != size (g)))
2832 {
2833 delete[] pEvalPoints;
2834 delete[] pMat;
2835 delete[] pL;
2836 delete[] coeffMonoms;
2837 delete[] pM;
2838 if (bufpEvalPoints != NULL)
2839 delete [] bufpEvalPoints;
2840 fail= true;
2841 if (alpha.level() != 1 && V_buf != alpha)
2842 prune1 (alpha);
2843 return 0;
2844 }
2845 CFIterator l= skel;
2846 for (CFIterator k= g; k.hasTerms(); k++, l++)
2847 {
2848 if (k.exp() != l.exp())
2849 {
2850 delete[] pEvalPoints;
2851 delete[] pMat;
2852 delete[] pL;
2853 delete[] coeffMonoms;
2854 delete[] pM;
2855 if (bufpEvalPoints != NULL)
2856 delete [] bufpEvalPoints;
2857 fail= true;
2858 if (alpha.level() != 1 && V_buf != alpha)
2859 prune1 (alpha);
2860 return 0;
2861 }
2862 }
2863 pEvalPoints[i + biggestSize]= evalPoints;
2864 gcds[i + biggestSize]= g;
2865 }
2866 }
2867 for (int i= 0; i < biggestSize; i++)
2868 {
2869 CFIterator l= gcds [i];
2870 evalPoints= pEvalPoints [i];
2871 for (int k= 1; k < skelSize; k++, l++)
2872 {
2873 if (i == 0)
2874 pMat[k]= CFMatrix (biggestSize2,coeffMonoms[k].size()+biggestSize2-1);
2875 if (k == minimalColumnsIndex)
2876 continue;
2877 coeffEval= evaluate (coeffMonoms[k], evalPoints);
2878 if (pMat[k].rows() >= i + 1)
2879 {
2880 for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2881 pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2882 }
2883 }
2884 }
2885 Mat= bufMat0;
2886 pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2 - 1);
2887 for (int j= 1; j <= Mat.rows(); j++)
2888 for (int k= 1; k <= Mat.columns(); k++)
2889 pMat [0] (j,k)= Mat (j,k);
2890 Mat= bufMat;
2891 for (int j= 1; j <= Mat.rows(); j++)
2892 for (int k= 1; k <= Mat.columns(); k++)
2893 pMat [minimalColumnsIndex] (j,k)= Mat (j,k);
2894 // write old matrix entries into new matrices
2895 for (int i= 0; i < skelSize; i++)
2896 {
2897 bufArray= pL[i];
2898 pL[i]= CFArray (biggestSize2);
2899 for (int j= 0; j < bufArray.size(); j++)
2900 pL[i] [j]= bufArray [j];
2901 }
2902 //write old vector entries into new and add new entries to old matrices
2903 for (int i= 0; i < biggestSize2 - biggestSize; i++)
2904 {
2905 CFIterator l= gcds [i + biggestSize];
2906 evalPoints= pEvalPoints [i + biggestSize];
2907 for (int k= 0; k < skelSize; k++, l++)
2908 {
2909 pL[k] [i + biggestSize]= l.coeff();
2910 coeffEval= evaluate (coeffMonoms[k], evalPoints);
2911 if (pMat[k].rows() >= i + biggestSize + 1)
2912 {
2913 for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2914 pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind - 1];
2915 }
2916 }
2917 }
2918 // begin new
2919 for (int i= 0; i < skelSize; i++)
2920 {
2921 if (pL[i].size() > 1)
2922 {
2923 for (int j= 2; j <= pMat[i].rows(); j++)
2924 pMat[i] (j, coeffMonoms[i].size() + j - 1)=
2925 -pL[i] [j - 1];
2926 }
2927 }
2928
2929 matColumns= biggestSize2 - 1;
2930 matRows= 0;
2931 for (int i= 0; i < skelSize; i++)
2932 {
2933 if (V_buf.level() == 1)
2934 (void) gaussianElimFp (pMat[i], pL[i]);
2935 else
2936 (void) gaussianElimFq (pMat[i], pL[i], V_buf);
2937
2938 if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)
2939 {
2940 delete[] pEvalPoints;
2941 delete[] pMat;
2942 delete[] pL;
2943 delete[] coeffMonoms;
2944 delete[] pM;
2945 if (bufpEvalPoints != NULL)
2946 delete [] bufpEvalPoints;
2947 fail= true;
2948 if (alpha.level() != 1 && V_buf != alpha)
2949 prune1 (alpha);
2950 return 0;
2951 }
2952 matRows += pMat[i].rows() - coeffMonoms[i].size();
2953 }
2954
2955 CFMatrix bufMat;
2956 Mat= CFMatrix (matRows, matColumns);
2957 ind= 0;
2958 bufArray= CFArray (matRows);
2959 CFArray bufArray2;
2960 for (int i= 0; i < skelSize; i++)
2961 {
2962 if (coeffMonoms[i].size() + 1 >= pMat[i].rows() || coeffMonoms[i].size() + 1 >= pMat[i].columns())
2963 {
2964 delete[] pEvalPoints;
2965 delete[] pMat;
2966 delete[] pL;
2967 delete[] coeffMonoms;
2968 delete[] pM;
2969 if (bufpEvalPoints != NULL)
2970 delete [] bufpEvalPoints;
2971 fail= true;
2972 if (alpha.level() != 1 && V_buf != alpha)
2973 prune1 (alpha);
2974 return 0;
2975 }
2976 bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),
2977 coeffMonoms[i].size() + 1, pMat[i].columns());
2978
2979 for (int j= 1; j <= bufMat.rows(); j++)
2980 for (int k= 1; k <= bufMat.columns(); k++)
2981 Mat (j + ind, k)= bufMat(j, k);
2982 bufArray2= coeffMonoms[i].size();
2983 for (int j= 1; j <= pMat[i].rows(); j++)
2984 {
2985 if (j > coeffMonoms[i].size())
2986 bufArray [j-coeffMonoms[i].size() + ind - 1]= pL[i] [j - 1];
2987 else
2988 bufArray2 [j - 1]= pL[i] [j - 1];
2989 }
2990 pL[i]= bufArray2;
2991 ind += bufMat.rows();
2992 pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());
2993 }
2994
2995 if (V_buf.level() != 1)
2996 solution= solveSystemFq (Mat, bufArray, V_buf);
2997 else
2998 solution= solveSystemFp (Mat, bufArray);
2999
3000 if (solution.size() == 0)
3001 {
3002 delete[] pEvalPoints;
3003 delete[] pMat;
3004 delete[] pL;
3005 delete[] coeffMonoms;
3006 delete[] pM;
3007 if (bufpEvalPoints != NULL)
3008 delete [] bufpEvalPoints;
3009 fail= true;
3010 if (alpha.level() != 1 && V_buf != alpha)
3011 prune1 (alpha);
3012 return 0;
3013 }
3014
3015 ind= 0;
3016 result= 0;
3017 CFArray bufSolution;
3018 for (int i= 0; i < skelSize; i++)
3019 {
3020 bufSolution= readOffSolution (pMat[i], pL[i], solution);
3021 for (int i= 0; i < bufSolution.size(); i++, ind++)
3022 result += Monoms [ind]*bufSolution[i];
3023 }
3024 if (alpha.level() != 1 && V_buf != alpha)
3025 {
3026 CFList u, v;
3027 result= mapDown (result,prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3028 prune1 (alpha);
3029 }
3030 result= N(result);
3031 delete[] pEvalPoints;
3032 delete[] pMat;
3033 delete[] pL;
3034 delete[] coeffMonoms;
3035 delete[] pM;
3036
3037 if (bufpEvalPoints != NULL)
3038 delete [] bufpEvalPoints;
3039 if (fdivides (result, F) && fdivides (result, G))
3040 return result;
3041 else
3042 {
3043 fail= true;
3044 return 0;
3045 }
3046 } // end of deKleine, Monagan & Wittkopf
3047
3048 result += Monoms[0];
3049 int ind2= 0, ind3= 2;
3050 ind= 0;
3051 for (int l= 0; l < minimalColumnsIndex; l++)
3052 ind += coeffMonoms[l].size();
3053 for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
3054 l++, ind2++, ind3++)
3055 {
3056 result += solution[l]*Monoms [1 + ind2];
3057 for (int i= 0; i < pMat[0].rows(); i++)
3058 firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
3059 }
3060 for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
3061 result += solution[l]*Monoms [ind + l];
3062 ind= coeffMonoms[0].size();
3063 for (int k= 1; k < skelSize; k++)
3064 {
3065 if (k == minimalColumnsIndex)
3066 {
3067 ind += coeffMonoms[k].size();
3068 continue;
3069 }
3070 if (k != minimalColumnsIndex)
3071 {
3072 for (int i= 0; i < biggestSize; i++)
3073 pL[k] [i] *= firstColumn [i];
3074 }
3075
3076 if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
3077 {
3078 bufArray= CFArray (coeffMonoms[k].size());
3079 for (int i= 0; i < bufArray.size(); i++)
3080 bufArray [i]= pL[k] [i];
3081 solution= solveGeneralVandermonde (pM[k], bufArray);
3082 }
3083 else
3084 solution= solveGeneralVandermonde (pM[k], pL[k]);
3085
3086 if (solution.size() == 0)
3087 {
3088 delete[] pEvalPoints;
3089 delete[] pMat;
3090 delete[] pL;
3091 delete[] coeffMonoms;
3092 delete[] pM;
3093 fail= true;
3094 if (alpha.level() != 1 && V_buf != alpha)
3095 prune1 (alpha);
3096 return 0;
3097 }
3098 if (k != minimalColumnsIndex)
3099 {
3100 for (int l= 0; l < solution.size(); l++)
3101 result += solution[l]*Monoms [ind + l];
3102 ind += solution.size();
3103 }
3104 }
3105
3106 delete[] pEvalPoints;
3107 delete[] pMat;
3108 delete[] pL;
3109 delete[] pM;
3110 delete[] coeffMonoms;
3111
3112 if (alpha.level() != 1 && V_buf != alpha)
3113 {
3114 CFList u, v;
3115 result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
3116 prune1 (alpha);
3117 }
3118 result= N(result);
3119
3120 if (fdivides (result, F) && fdivides (result, G))
3121 return result;
3122 else
3123 {
3124 fail= true;
3125 return 0;
3126 }
3127}
3128
3130 const Variable & alpha, CFList& l, bool& topLevel)
3131{
3132 CanonicalForm A= F;
3133 CanonicalForm B= G;
3134 if (F.isZero() && degree(G) > 0) return G/Lc(G);
3135 else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3136 else if (F.isZero() && G.isZero()) return F.genOne();
3137 if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3138 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3139 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3140 if (F == G) return F/Lc(F);
3141
3142 CFMap M,N;
3143 int best_level= myCompress (A, B, M, N, topLevel);
3144
3145 if (best_level == 0) return B.genOne();
3146
3147 A= M(A);
3148 B= M(B);
3149
3150 Variable x= Variable (1);
3151
3152 //univariate case
3153 if (A.isUnivariate() && B.isUnivariate())
3154 return N (gcd (A, B));
3155
3156 CanonicalForm cA, cB; // content of A and B
3157 CanonicalForm ppA, ppB; // primitive part of A and B
3158 CanonicalForm gcdcAcB;
3159
3160 cA = uni_content (A);
3161 cB = uni_content (B);
3162 gcdcAcB= gcd (cA, cB);
3163 ppA= A/cA;
3164 ppB= B/cB;
3165
3166 CanonicalForm lcA, lcB; // leading coefficients of A and B
3167 CanonicalForm gcdlcAlcB;
3168 lcA= uni_lcoeff (ppA);
3169 lcB= uni_lcoeff (ppB);
3170
3171 if (fdivides (lcA, lcB))
3172 {
3173 if (fdivides (A, B))
3174 return F/Lc(F);
3175 }
3176 if (fdivides (lcB, lcA))
3177 {
3178 if (fdivides (B, A))
3179 return G/Lc(G);
3180 }
3181
3182 gcdlcAlcB= gcd (lcA, lcB);
3183
3184 int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3185 int d0;
3186
3187 if (d == 0)
3188 return N(gcdcAcB);
3189 d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3190
3191 if (d0 < d)
3192 d= d0;
3193
3194 if (d == 0)
3195 return N(gcdcAcB);
3196
3197 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3198 CanonicalForm newtonPoly= 1;
3199 m= gcdlcAlcB;
3200 G_m= 0;
3201 H= 0;
3202 bool fail= false;
3203 topLevel= false;
3204 bool inextension= false;
3205 Variable V_buf= alpha, V_buf4= alpha;
3206 CanonicalForm prim_elem, im_prim_elem;
3207 CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
3208 CFList source, dest;
3209 do // first do
3210 {
3211 random_element= randomElement (m, V_buf, l, fail);
3212 if (random_element == 0 && !fail)
3213 {
3214 if (!find (l, random_element))
3215 l.append (random_element);
3216 continue;
3217 }
3218 if (fail)
3219 {
3220 source= CFList();
3221 dest= CFList();
3222
3223 Variable V_buf3= V_buf;
3224 V_buf= chooseExtension (V_buf);
3225 bool prim_fail= false;
3226 Variable V_buf2;
3227 prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
3228 if (V_buf4 == alpha)
3229 prim_elem_alpha= prim_elem;
3230
3231 if (V_buf3 != V_buf4)
3232 {
3233 m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3234 G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3235 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3236 source, dest);
3237 ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3238 ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3239 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4, source,
3240 dest);
3241 for (CFListIterator i= l; i.hasItem(); i++)
3242 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3243 source, dest);
3244 }
3245
3246 ASSERT (!prim_fail, "failure in integer factorizer");
3247 if (prim_fail)
3248 ; //ERROR
3249 else
3250 im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
3251
3252 if (V_buf4 == alpha)
3253 im_prim_elem_alpha= im_prim_elem;
3254 else
3255 im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
3256 im_prim_elem, source, dest);
3257
3258 DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3259 DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3260 inextension= true;
3261 for (CFListIterator i= l; i.hasItem(); i++)
3262 i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3263 im_prim_elem, source, dest);
3264 m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3265 G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3266 newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3267 source, dest);
3268 ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3269 ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3270 gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3271 source, dest);
3272
3273 fail= false;
3274 random_element= randomElement (m, V_buf, l, fail );
3275 DEBOUTLN (cerr, "fail= " << fail);
3276 CFList list;
3277 TIMING_START (gcd_recursion);
3278 G_random_element=
3279 sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3280 list, topLevel);
3281 TIMING_END_AND_PRINT (gcd_recursion,
3282 "time for recursive call: ");
3283 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3284 V_buf4= V_buf;
3285 }
3286 else
3287 {
3288 CFList list;
3289 TIMING_START (gcd_recursion);
3290 G_random_element=
3291 sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
3292 list, topLevel);
3293 TIMING_END_AND_PRINT (gcd_recursion,
3294 "time for recursive call: ");
3295 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3296 }
3297
3298 if (!G_random_element.inCoeffDomain())
3299 d0= totaldegree (G_random_element, Variable(2),
3300 Variable (G_random_element.level()));
3301 else
3302 d0= 0;
3303
3304 if (d0 == 0)
3305 {
3306 if (inextension)
3307 prune1 (alpha);
3308 return N(gcdcAcB);
3309 }
3310 if (d0 > d)
3311 {
3312 if (!find (l, random_element))
3313 l.append (random_element);
3314 continue;
3315 }
3316
3317 G_random_element=
3318 (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3319 * G_random_element;
3320
3321 skeleton= G_random_element;
3322 if (!G_random_element.inCoeffDomain())
3323 d0= totaldegree (G_random_element, Variable(2),
3324 Variable (G_random_element.level()));
3325 else
3326 d0= 0;
3327
3328 if (d0 < d)
3329 {
3330 m= gcdlcAlcB;
3331 newtonPoly= 1;
3332 G_m= 0;
3333 d= d0;
3334 }
3335
3336 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3337 if (uni_lcoeff (H) == gcdlcAlcB)
3338 {
3339 cH= uni_content (H);
3340 ppH= H/cH;
3341 if (inextension)
3342 {
3343 CFList u, v;
3344 //maybe it's better to test if ppH is an element of F(\alpha) before
3345 //mapping down
3346 if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3347 {
3348 DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3349 ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3350 ppH /= Lc(ppH);
3351 DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3352 prune1 (alpha);
3353 return N(gcdcAcB*ppH);
3354 }
3355 }
3356 else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3357 return N(gcdcAcB*ppH);
3358 }
3359 G_m= H;
3360 newtonPoly= newtonPoly*(x - random_element);
3361 m= m*(x - random_element);
3362 if (!find (l, random_element))
3363 l.append (random_element);
3364
3365 if (getCharacteristic () > 3 && size (skeleton) < 100)
3366 {
3367 CFArray Monoms;
3368 CFArray *coeffMonoms;
3369 do //second do
3370 {
3371 random_element= randomElement (m, V_buf, l, fail);
3372 if (random_element == 0 && !fail)
3373 {
3374 if (!find (l, random_element))
3375 l.append (random_element);
3376 continue;
3377 }
3378 if (fail)
3379 {
3380 source= CFList();
3381 dest= CFList();
3382
3383 Variable V_buf3= V_buf;
3384 V_buf= chooseExtension (V_buf);
3385 bool prim_fail= false;
3386 Variable V_buf2;
3387 prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
3388 if (V_buf4 == alpha)
3389 prim_elem_alpha= prim_elem;
3390
3391 if (V_buf3 != V_buf4)
3392 {
3393 m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3394 G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3395 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3396 source, dest);
3397 ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3398 ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3399 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
3400 source, dest);
3401 for (CFListIterator i= l; i.hasItem(); i++)
3402 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3403 source, dest);
3404 }
3405
3406 ASSERT (!prim_fail, "failure in integer factorizer");
3407 if (prim_fail)
3408 ; //ERROR
3409 else
3410 im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
3411
3412 if (V_buf4 == alpha)
3413 im_prim_elem_alpha= im_prim_elem;
3414 else
3415 im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
3416 prim_elem, im_prim_elem, source, dest);
3417
3418 DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3419 DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3420 inextension= true;
3421 for (CFListIterator i= l; i.hasItem(); i++)
3422 i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3423 im_prim_elem, source, dest);
3424 m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3425 G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3426 newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3427 source, dest);
3428 ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3429 ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3430
3431 gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3432 source, dest);
3433
3434 fail= false;
3435 random_element= randomElement (m, V_buf, l, fail);
3436 DEBOUTLN (cerr, "fail= " << fail);
3437 CFList list;
3438 TIMING_START (gcd_recursion);
3439
3440 V_buf4= V_buf;
3441
3442 //sparseInterpolation
3443 bool sparseFail= false;
3444 if (LC (skeleton).inCoeffDomain())
3445 G_random_element=
3446 monicSparseInterpol (ppA (random_element, x), ppB(random_element,x),
3447 skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
3448 else
3449 G_random_element=
3450 nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3451 skeleton, V_buf, sparseFail, coeffMonoms,
3452 Monoms);
3453 if (sparseFail)
3454 break;
3455
3456 TIMING_END_AND_PRINT (gcd_recursion,
3457 "time for recursive call: ");
3458 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3459 }
3460 else
3461 {
3462 CFList list;
3463 TIMING_START (gcd_recursion);
3464 bool sparseFail= false;
3465 if (LC (skeleton).inCoeffDomain())
3466 G_random_element=
3467 monicSparseInterpol (ppA (random_element, x),ppB(random_element, x),
3468 skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
3469 else
3470 G_random_element=
3471 nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3472 skeleton, V_buf, sparseFail, coeffMonoms,
3473 Monoms);
3474 if (sparseFail)
3475 break;
3476
3477 TIMING_END_AND_PRINT (gcd_recursion,
3478 "time for recursive call: ");
3479 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3480 }
3481
3482 if (!G_random_element.inCoeffDomain())
3483 d0= totaldegree (G_random_element, Variable(2),
3484 Variable (G_random_element.level()));
3485 else
3486 d0= 0;
3487
3488 if (d0 == 0)
3489 {
3490 if (inextension)
3491 prune1 (alpha);
3492 return N(gcdcAcB);
3493 }
3494 if (d0 > d)
3495 {
3496 if (!find (l, random_element))
3497 l.append (random_element);
3498 continue;
3499 }
3500
3501 G_random_element=
3502 (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3503 * G_random_element;
3504
3505 if (!G_random_element.inCoeffDomain())
3506 d0= totaldegree (G_random_element, Variable(2),
3507 Variable (G_random_element.level()));
3508 else
3509 d0= 0;
3510
3511 if (d0 < d)
3512 {
3513 m= gcdlcAlcB;
3514 newtonPoly= 1;
3515 G_m= 0;
3516 d= d0;
3517 }
3518
3519 TIMING_START (newton_interpolation);
3520 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3521 TIMING_END_AND_PRINT (newton_interpolation,
3522 "time for newton interpolation: ");
3523
3524 //termination test
3525 if (uni_lcoeff (H) == gcdlcAlcB)
3526 {
3527 cH= uni_content (H);
3528 ppH= H/cH;
3529 if (inextension)
3530 {
3531 CFList u, v;
3532 //maybe it's better to test if ppH is an element of F(\alpha) before
3533 //mapping down
3534 if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3535 {
3536 DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3537 ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3538 ppH /= Lc(ppH);
3539 DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3540 prune1 (alpha);
3541 return N(gcdcAcB*ppH);
3542 }
3543 }
3544 else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3545 {
3546 return N(gcdcAcB*ppH);
3547 }
3548 }
3549
3550 G_m= H;
3551 newtonPoly= newtonPoly*(x - random_element);
3552 m= m*(x - random_element);
3553 if (!find (l, random_element))
3554 l.append (random_element);
3555
3556 } while (1);
3557 }
3558 } while (1);
3559}
3560
3562 bool& topLevel, CFList& l)
3563{
3564 CanonicalForm A= F;
3565 CanonicalForm B= G;
3566 if (F.isZero() && degree(G) > 0) return G/Lc(G);
3567 else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3568 else if (F.isZero() && G.isZero()) return F.genOne();
3569 if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3570 if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3571 if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3572 if (F == G) return F/Lc(F);
3573
3574 CFMap M,N;
3575 int best_level= myCompress (A, B, M, N, topLevel);
3576
3577 if (best_level == 0) return B.genOne();
3578
3579 A= M(A);
3580 B= M(B);
3581
3582 Variable x= Variable (1);
3583
3584 //univariate case
3585 if (A.isUnivariate() && B.isUnivariate())
3586 return N (gcd (A, B));
3587
3588 CanonicalForm cA, cB; // content of A and B
3589 CanonicalForm ppA, ppB; // primitive part of A and B
3590 CanonicalForm gcdcAcB;
3591
3592 cA = uni_content (A);
3593 cB = uni_content (B);
3594 gcdcAcB= gcd (cA, cB);
3595 ppA= A/cA;
3596 ppB= B/cB;
3597
3598 CanonicalForm lcA, lcB; // leading coefficients of A and B
3599 CanonicalForm gcdlcAlcB;
3600 lcA= uni_lcoeff (ppA);
3601 lcB= uni_lcoeff (ppB);
3602
3603 if (fdivides (lcA, lcB))
3604 {
3605 if (fdivides (A, B))
3606 return F/Lc(F);
3607 }
3608 if (fdivides (lcB, lcA))
3609 {
3610 if (fdivides (B, A))
3611 return G/Lc(G);
3612 }
3613
3614 gcdlcAlcB= gcd (lcA, lcB);
3615
3616 int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3617 int d0;
3618
3619 if (d == 0)
3620 return N(gcdcAcB);
3621
3622 d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3623
3624 if (d0 < d)
3625 d= d0;
3626
3627 if (d == 0)
3628 return N(gcdcAcB);
3629
3630 CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3631 CanonicalForm newtonPoly= 1;
3632 m= gcdlcAlcB;
3633 G_m= 0;
3634 H= 0;
3635 bool fail= false;
3636 topLevel= false;
3637 bool inextension= false;
3638 Variable V_buf, alpha, cleanUp;
3639 CanonicalForm prim_elem, im_prim_elem;
3640 CFList source, dest;
3641 do //first do
3642 {
3643 if (inextension)
3644 random_element= randomElement (m, V_buf, l, fail);
3645 else
3646 random_element= FpRandomElement (m, l, fail);
3647 if (random_element == 0 && !fail)
3648 {
3649 if (!find (l, random_element))
3650 l.append (random_element);
3651 continue;
3652 }
3653
3654 if (!fail && !inextension)
3655 {
3656 CFList list;
3657 TIMING_START (gcd_recursion);
3658 G_random_element=
3659 sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
3660 list);
3661 TIMING_END_AND_PRINT (gcd_recursion,
3662 "time for recursive call: ");
3663 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3664 }
3665 else if (!fail && inextension)
3666 {
3667 CFList list;
3668 TIMING_START (gcd_recursion);
3669 G_random_element=
3670 sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3671 list, topLevel);
3672 TIMING_END_AND_PRINT (gcd_recursion,
3673 "time for recursive call: ");
3674 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3675 }
3676 else if (fail && !inextension)
3677 {
3678 source= CFList();
3679 dest= CFList();
3680 CFList list;
3682 int deg= 2;
3683 bool initialized= false;
3684 do
3685 {
3686 mipo= randomIrredpoly (deg, x);
3687 if (initialized)
3688 setMipo (alpha, mipo);
3689 else
3690 alpha= rootOf (mipo);
3691 inextension= true;
3692 fail= false;
3693 initialized= true;
3694 random_element= randomElement (m, alpha, l, fail);
3695 deg++;
3696 } while (fail);
3697 cleanUp= alpha;
3698 V_buf= alpha;
3699 list= CFList();
3700 TIMING_START (gcd_recursion);
3701 G_random_element=
3702 sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3703 list, topLevel);
3704 TIMING_END_AND_PRINT (gcd_recursion,
3705 "time for recursive call: ");
3706 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3707 }
3708 else if (fail && inextension)
3709 {
3710 source= CFList();
3711 dest= CFList();
3712
3713 Variable V_buf3= V_buf;
3714 V_buf= chooseExtension (V_buf);
3715 bool prim_fail= false;
3716 Variable V_buf2;
3717 CanonicalForm prim_elem, im_prim_elem;
3718 prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3719
3720 if (V_buf3 != alpha)
3721 {
3722 m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3723 G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3724 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source,
3725 dest);
3726 ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3727 ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3728 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3729 dest);
3730 for (CFListIterator i= l; i.hasItem(); i++)
3731 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3732 source, dest);
3733 }
3734
3735 ASSERT (!prim_fail, "failure in integer factorizer");
3736 if (prim_fail)
3737 ; //ERROR
3738 else
3739 im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3740
3741 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3742 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3743
3744 for (CFListIterator i= l; i.hasItem(); i++)
3745 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3746 im_prim_elem, source, dest);
3747 m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3748 G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3749 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3750 source, dest);
3751 ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3752 ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3753 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3754 source, dest);
3755 fail= false;
3756 random_element= randomElement (m, V_buf, l, fail );
3757 DEBOUTLN (cerr, "fail= " << fail);
3758 CFList list;
3759 TIMING_START (gcd_recursion);
3760 G_random_element=
3761 sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3762 list, topLevel);
3763 TIMING_END_AND_PRINT (gcd_recursion,
3764 "time for recursive call: ");
3765 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3766 alpha= V_buf;
3767 }
3768
3769 if (!G_random_element.inCoeffDomain())
3770 d0= totaldegree (G_random_element, Variable(2),
3771 Variable (G_random_element.level()));
3772 else
3773 d0= 0;
3774
3775 if (d0 == 0)
3776 {
3777 if (inextension)
3778 prune (cleanUp);
3779 return N(gcdcAcB);
3780 }
3781 if (d0 > d)
3782 {
3783 if (!find (l, random_element))
3784 l.append (random_element);
3785 continue;
3786 }
3787
3788 G_random_element=
3789 (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3790 * G_random_element;
3791
3792 skeleton= G_random_element;
3793
3794 if (!G_random_element.inCoeffDomain())
3795 d0= totaldegree (G_random_element, Variable(2),
3796 Variable (G_random_element.level()));
3797 else
3798 d0= 0;
3799
3800 if (d0 < d)
3801 {
3802 m= gcdlcAlcB;
3803 newtonPoly= 1;
3804 G_m= 0;
3805 d= d0;
3806 }
3807
3808 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3809
3810 if (uni_lcoeff (H) == gcdlcAlcB)
3811 {
3812 cH= uni_content (H);
3813 ppH= H/cH;
3814 ppH /= Lc (ppH);
3815 DEBOUTLN (cerr, "ppH= " << ppH);
3816
3817 if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3818 {
3819 if (inextension)
3820 prune (cleanUp);
3821 return N(gcdcAcB*ppH);
3822 }
3823 }
3824 G_m= H;
3825 newtonPoly= newtonPoly*(x - random_element);
3826 m= m*(x - random_element);
3827 if (!find (l, random_element))
3828 l.append (random_element);
3829
3830 if ((getCharacteristic() > 3 && size (skeleton) < 200))
3831 {
3832 CFArray Monoms;
3833 CFArray* coeffMonoms;
3834
3835 do //second do
3836 {
3837 if (inextension)
3838 random_element= randomElement (m, alpha, l, fail);
3839 else
3840 random_element= FpRandomElement (m, l, fail);
3841 if (random_element == 0 && !fail)
3842 {
3843 if (!find (l, random_element))
3844 l.append (random_element);
3845 continue;
3846 }
3847
3848 bool sparseFail= false;
3849 if (!fail && !inextension)
3850 {
3851 CFList list;
3852 TIMING_START (gcd_recursion);
3853 if (LC (skeleton).inCoeffDomain())
3854 G_random_element=
3855 monicSparseInterpol(ppA(random_element, x), ppB (random_element, x),
3856 skeleton, x, sparseFail, coeffMonoms,
3857 Monoms);
3858 else
3859 G_random_element=
3860 nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3861 skeleton, x, sparseFail,
3862 coeffMonoms, Monoms);
3863 TIMING_END_AND_PRINT (gcd_recursion,
3864 "time for recursive call: ");
3865 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3866 }
3867 else if (!fail && inextension)
3868 {
3869 CFList list;
3870 TIMING_START (gcd_recursion);
3871 if (LC (skeleton).inCoeffDomain())
3872 G_random_element=
3873 monicSparseInterpol(ppA (random_element,x), ppB (random_element, x),
3874 skeleton, alpha, sparseFail, coeffMonoms,
3875 Monoms);
3876 else
3877 G_random_element=
3878 nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3879 skeleton, alpha, sparseFail, coeffMonoms,
3880 Monoms);
3881 TIMING_END_AND_PRINT (gcd_recursion,
3882 "time for recursive call: ");
3883 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3884 }
3885 else if (fail && !inextension)
3886 {
3887 source= CFList();
3888 dest= CFList();
3889 CFList list;
3891 int deg= 2;
3892 bool initialized= false;
3893 do
3894 {
3895 mipo= randomIrredpoly (deg, x);
3896 if (initialized)
3897 setMipo (alpha, mipo);
3898 else
3899 alpha= rootOf (mipo);
3900 inextension= true;
3901 fail= false;
3902 initialized= true;
3903 random_element= randomElement (m, alpha, l, fail);
3904 deg++;
3905 } while (fail);
3906 cleanUp= alpha;
3907 V_buf= alpha;
3908 list= CFList();
3909 TIMING_START (gcd_recursion);
3910 if (LC (skeleton).inCoeffDomain())
3911 G_random_element=
3912 monicSparseInterpol (ppA (random_element,x), ppB (random_element,x),
3913 skeleton, alpha, sparseFail, coeffMonoms,
3914 Monoms);
3915 else
3916 G_random_element=
3917 nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3918 skeleton, alpha, sparseFail, coeffMonoms,
3919 Monoms);
3920 TIMING_END_AND_PRINT (gcd_recursion,
3921 "time for recursive call: ");
3922 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3923 }
3924 else if (fail && inextension)
3925 {
3926 source= CFList();
3927 dest= CFList();
3928
3929 Variable V_buf3= V_buf;
3930 V_buf= chooseExtension (V_buf);
3931 bool prim_fail= false;
3932 Variable V_buf2;
3933 CanonicalForm prim_elem, im_prim_elem;
3934 prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3935
3936 if (V_buf3 != alpha)
3937 {
3938 m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3939 G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3940 newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3941 source, dest);
3942 ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3943 ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3944 gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3945 source, dest);
3946 for (CFListIterator i= l; i.hasItem(); i++)
3947 i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3948 source, dest);
3949 }
3950
3951 ASSERT (!prim_fail, "failure in integer factorizer");
3952 if (prim_fail)
3953 ; //ERROR
3954 else
3955 im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3956
3957 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3958 DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3959
3960 for (CFListIterator i= l; i.hasItem(); i++)
3961 i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3962 im_prim_elem, source, dest);
3963 m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3964 G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3965 newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3966 source, dest);
3967 ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3968 ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3969 gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3970 source, dest);
3971 fail= false;
3972 random_element= randomElement (m, V_buf, l, fail );
3973 DEBOUTLN (cerr, "fail= " << fail);
3974 CFList list;
3975 TIMING_START (gcd_recursion);
3976 if (LC (skeleton).inCoeffDomain())
3977 G_random_element=
3978 monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
3979 skeleton, V_buf, sparseFail, coeffMonoms,
3980 Monoms);
3981 else
3982 G_random_element=
3983 nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
3984 skeleton, V_buf, sparseFail,
3985 coeffMonoms, Monoms);
3986 TIMING_END_AND_PRINT (gcd_recursion,
3987 "time for recursive call: ");
3988 DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3989 alpha= V_buf;
3990 }
3991
3992 if (sparseFail)
3993 break;
3994
3995 if (!G_random_element.inCoeffDomain())
3996 d0= totaldegree (G_random_element, Variable(2),
3997 Variable (G_random_element.level()));
3998 else
3999 d0= 0;
4000
4001 if (d0 == 0)
4002 {
4003 if (inextension)
4004 prune (cleanUp);
4005 return N(gcdcAcB);
4006 }
4007 if (d0 > d)
4008 {
4009 if (!find (l, random_element))
4010 l.append (random_element);
4011 continue;
4012 }
4013
4014 G_random_element=
4015 (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
4016 * G_random_element;
4017
4018 if (!G_random_element.inCoeffDomain())
4019 d0= totaldegree (G_random_element, Variable(2),
4020 Variable (G_random_element.level()));
4021 else
4022 d0= 0;
4023
4024 if (d0 < d)
4025 {
4026 m= gcdlcAlcB;
4027 newtonPoly= 1;
4028 G_m= 0;
4029 d= d0;
4030 }
4031
4032 TIMING_START (newton_interpolation);
4033 H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
4034 TIMING_END_AND_PRINT (newton_interpolation,
4035 "time for newton interpolation: ");
4036
4037 //termination test
4038 if (uni_lcoeff (H) == gcdlcAlcB)
4039 {
4040 cH= uni_content (H);
4041 ppH= H/cH;
4042 ppH /= Lc (ppH);
4043 DEBOUTLN (cerr, "ppH= " << ppH);
4044 if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
4045 {
4046 if (inextension)
4047 prune (cleanUp);
4048 return N(gcdcAcB*ppH);
4049 }
4050 }
4051
4052 G_m= H;
4053 newtonPoly= newtonPoly*(x - random_element);
4054 m= m*(x - random_element);
4055 if (!find (l, random_element))
4056 l.append (random_element);
4057
4058 } while (1); //end of second do
4059 }
4060 else
4061 {
4062 if (inextension)
4063 prune (cleanUp);
4064 return N(gcdcAcB*modGCDFp (ppA, ppB));
4065 }
4066 } while (1); //end of first do
4067}
4068
4069TIMING_DEFINE_PRINT(modZ_termination)
4070TIMING_DEFINE_PRINT(modZ_recursion)
4071
4072/// modular gcd algorithm, see Keith, Czapor, Geddes "Algorithms for Computer
4073/// Algebra", Algorithm 7.1
4075{
4076 CanonicalForm f, g, cl, q(0), Dp, newD, D, newq, newqh;
4077 int p, i, dp_deg, d_deg=-1;
4078
4079 CanonicalForm cd ( bCommonDen( FF ));
4080 f=cd*FF;
4083 cf= icontent (f);
4084 f /= cf;
4085 //cd = bCommonDen( f ); f *=cd;
4086 //f /=vcontent(f,Variable(1));
4087
4090 cg= icontent (g);
4091 g /= cg;
4092 //cd = bCommonDen( g ); g *=cd;
4093 //g /=vcontent(g,Variable(1));
4094
4097 lcf= f.lc();
4098 lcg= g.lc();
4099 cl = gcd (f.lc(),g.lc());
4104 for (i= tmax (f.level(), g.level()); i > 0; i--)
4105 {
4106 if (degree (f, i) <= 0 || degree (g, i) <= 0)
4107 continue;
4108 else
4109 {
4110 minCommonDeg= tmin (degree (g, i), degree (f, i));
4111 break;
4112 }
4113 }
4114 if (i == 0)
4115 return gcdcfcg;
4116 for (; i > 0; i--)
4117 {
4118 if (degree (f, i) <= 0 || degree (g, i) <= 0)
4119 continue;
4120 else
4122 }
4123 b= 2*tmin (maxNorm (f), maxNorm (g))*abs (cl)*
4125 bool equal= false;
4126 i = cf_getNumBigPrimes() - 1;
4127
4130 //Off (SW_RATIONAL);
4131 while ( true )
4132 {
4133 p = cf_getBigPrime( i );
4134 i--;
4135 while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 )
4136 {
4137 p = cf_getBigPrime( i );
4138 i--;
4139 }
4140 //printf("try p=%d\n",p);
4142 fp= mapinto (f);
4143 gp= mapinto (g);
4144 TIMING_START (modZ_recursion)
4145#if defined(HAVE_NTL) || defined(HAVE_FLINT)
4146 if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500)
4147 Dp = modGCDFp (fp, gp, cofp, cogp);
4148 else
4149 {
4150 Dp= gcd_poly (fp, gp);
4151 cofp= fp/Dp;
4152 cogp= gp/Dp;
4153 }
4154#else
4155 Dp= gcd_poly (fp, gp);
4156 cofp= fp/Dp;
4157 cogp= gp/Dp;
4158#endif
4159 TIMING_END_AND_PRINT (modZ_recursion,
4160 "time for gcd mod p in modular gcd: ");
4161 Dp /=Dp.lc();
4162 Dp *= mapinto (cl);
4163 cofp /= lc (cofp);
4164 cofp *= mapinto (lcf);
4165 cogp /= lc (cogp);
4166 cogp *= mapinto (lcg);
4167 setCharacteristic( 0 );
4168 dp_deg=totaldegree(Dp);
4169 if ( dp_deg == 0 )
4170 {
4171 //printf(" -> 1\n");
4172 return CanonicalForm(gcdcfcg);
4173 }
4174 if ( q.isZero() )
4175 {
4176 D = mapinto( Dp );
4177 cof= mapinto (cofp);
4178 cog= mapinto (cogp);
4179 d_deg=dp_deg;
4180 q = p;
4181 Dn= balance_p (D, p);
4182 cofn= balance_p (cof, p);
4183 cogn= balance_p (cog, p);
4184 }
4185 else
4186 {
4187 if ( dp_deg == d_deg )
4188 {
4189 chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
4190 chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
4191 chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
4192 cof= newCof;
4193 cog= newCog;
4194 newqh= newq/2;
4195 Dn= balance_p (newD, newq, newqh);
4196 cofn= balance_p (newCof, newq, newqh);
4197 cogn= balance_p (newCog, newq, newqh);
4198 if (test != Dn) //balance_p (newD, newq))
4199 test= balance_p (newD, newq);
4200 else
4201 equal= true;
4202 q = newq;
4203 D = newD;
4204 }
4205 else if ( dp_deg < d_deg )
4206 {
4207 //printf(" were all bad, try more\n");
4208 // all previous p's are bad primes
4209 q = p;
4210 D = mapinto( Dp );
4211 cof= mapinto (cof);
4212 cog= mapinto (cog);
4213 d_deg=dp_deg;
4214 test= 0;
4215 equal= false;
4216 Dn= balance_p (D, p);
4217 cofn= balance_p (cof, p);
4218 cogn= balance_p (cog, p);
4219 }
4220 else
4221 {
4222 //printf(" was bad, try more\n");
4223 }
4224 //else dp_deg > d_deg: bad prime
4225 }
4226 if ( i >= 0 )
4227 {
4228 cDn= icontent (Dn);
4229 Dn /= cDn;
4230 cofn /= cl/cDn;
4231 //cofn /= icontent (cofn);
4232 cogn /= cl/cDn;
4233 //cogn /= icontent (cogn);
4234 TIMING_START (modZ_termination);
4235 if ((terminationTest (f,g, cofn, cogn, Dn)) ||
4236 ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
4237 {
4238 TIMING_END_AND_PRINT (modZ_termination,
4239 "time for successful termination in modular gcd: ");
4240 //printf(" -> success\n");
4241 return Dn*gcdcfcg;
4242 }
4243 TIMING_END_AND_PRINT (modZ_termination,
4244 "time for unsuccessful termination in modular gcd: ");
4245 equal= false;
4246 //else: try more primes
4247 }
4248 else
4249 { // try other method
4250 //printf("try other gcd\n");
4252 D=gcd_poly( f, g );
4254 return D*gcdcfcg;
4255 }
4256 }
4257}
4258#endif
CFMatrix * convertNmod_mat_t2FacCFMatrix(const nmod_mat_t m)
conversion of a FLINT matrix over Z/p to a factory matrix
void convertFacCFMatrix2nmod_mat_t(nmod_mat_t M, const CFMatrix &m)
conversion of a factory matrix over Z/p to a nmod_mat_t
void convertFacCFMatrix2Fq_nmod_mat_t(fq_nmod_mat_t M, const fq_nmod_ctx_t fq_con, const CFMatrix &m)
conversion of a factory matrix over F_q to a fq_nmod_mat_t
CanonicalForm convertnmod_poly_t2FacCF(const nmod_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z/p to CanonicalForm
CFMatrix * convertFq_nmod_mat_t2FacCFMatrix(const fq_nmod_mat_t m, const fq_nmod_ctx_t &fq_con, const Variable &alpha)
conversion of a FLINT matrix over F_q to a factory matrix
This file defines functions for conversion to FLINT (www.flintlib.org) and back.
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:411
Rational abs(const Rational &a)
Definition: GMPrat.cc:436
CFMatrix * convertNTLmat_zz_p2FacCFMatrix(const mat_zz_p &m)
Definition: NTLconvert.cc:1183
CFMatrix * convertNTLmat_zz_pE2FacCFMatrix(const mat_zz_pE &m, const Variable &alpha)
Definition: NTLconvert.cc:1212
CanonicalForm convertNTLzzpX2CF(const zz_pX &poly, const Variable &x)
Definition: NTLconvert.cc:255
mat_zz_pE * convertFacCFMatrix2NTLmat_zz_pE(const CFMatrix &m)
Definition: NTLconvert.cc:1196
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
Definition: NTLconvert.cc:105
mat_zz_p * convertFacCFMatrix2NTLmat_zz_p(const CFMatrix &m)
Definition: NTLconvert.cc:1167
VAR long fac_NTL_char
Definition: NTLconvert.cc:46
Conversion to and from NTL.
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 FACTORY_PUBLIC content(const CanonicalForm &)
CanonicalForm content ( const CanonicalForm & f )
Definition: cf_gcd.cc:603
CanonicalForm mapinto(const CanonicalForm &f)
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
int getNumVars(const CanonicalForm &f)
int getNumVars ( const CanonicalForm & f )
Definition: cf_ops.cc:314
CanonicalForm lc(const CanonicalForm &f)
CanonicalForm FACTORY_PUBLIC icontent(const CanonicalForm &f)
CanonicalForm icontent ( const CanonicalForm & f )
Definition: cf_gcd.cc:74
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)
int getGFDegree()
Definition: cf_char.cc:75
Array< CanonicalForm > CFArray
void FACTORY_PUBLIC setCharacteristic(int c)
Definition: cf_char.cc:28
Matrix< CanonicalForm > CFMatrix
CanonicalForm FACTORY_PUBLIC swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition: cf_ops.cc:168
int totaldegree(const CanonicalForm &f)
int totaldegree ( const CanonicalForm & f )
Definition: cf_ops.cc:523
CanonicalForm Lc(const CanonicalForm &f)
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
List< CanonicalForm > CFList
CanonicalForm LC(const CanonicalForm &f)
int FACTORY_PUBLIC getCharacteristic()
Definition: cf_char.cc:70
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:56
int * degsf
Definition: cfEzgcd.cc:59
int f_zero
Definition: cfEzgcd.cc:69
int l
Definition: cfEzgcd.cc:100
int m
Definition: cfEzgcd.cc:128
const CanonicalForm CFMap CFMap int & both_non_zero
Definition: cfEzgcd.cc:57
int g_zero
Definition: cfEzgcd.cc:70
int k
Definition: cfEzgcd.cc:99
int * degsg
Definition: cfEzgcd.cc:60
const CanonicalForm CFMap CFMap bool topLevel
Definition: cfGcdAlgExt.cc:57
int both_zero
Definition: cfGcdAlgExt.cc:71
coprimality check and change of representation mod n
CanonicalForm sparseGCDFq(const CanonicalForm &F, const CanonicalForm &G, const Variable &alpha, CFList &l, bool &topLevel)
Definition: cfModGcd.cc:3129
CFArray readOffSolution(const CFMatrix &M, const long rk)
M in row echolon form, rk rank of M.
Definition: cfModGcd.cc:1688
CanonicalForm modGCDFq(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, Variable &alpha, CFList &l, bool &topLevel)
GCD of F and G over , l and topLevel are only used internally, output is monic based on Alg....
Definition: cfModGcd.cc:478
static CanonicalForm GFRandomElement(const CanonicalForm &F, CFList &list, bool &fail)
compute a random element a of GF, s.t. F(a) , F is a univariate polynomial, returns fail if there ar...
Definition: cfModGcd.cc:819
CFArray evaluateMonom(const CanonicalForm &F, const CFList &evalPoints)
Definition: cfModGcd.cc:1991
const CanonicalForm const CanonicalForm const CanonicalForm & coG
Definition: cfModGcd.cc:67
CFArray getMonoms(const CanonicalForm &F)
extract monomials of F, parts in algebraic variable are considered coefficients
Definition: cfModGcd.cc:1956
void mult(CFList &L1, const CFList &L2)
multiply two lists componentwise
Definition: cfModGcd.cc:2175
const CanonicalForm const CanonicalForm const CanonicalForm const CanonicalForm & cand
Definition: cfModGcd.cc:69
CanonicalForm cofp
Definition: cfModGcd.cc:4128
long gaussianElimFq(CFMatrix &M, CFArray &L, const Variable &alpha)
Definition: cfModGcd.cc:1784
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:4081
CanonicalForm lcg
Definition: cfModGcd.cc:4096
static Variable chooseExtension(const Variable &alpha)
Definition: cfModGcd.cc:420
int dp_deg
Definition: cfModGcd.cc:4077
CanonicalForm newCog
Definition: cfModGcd.cc:4128
CFArray evaluate(const CFArray &A, const CFList &evalPoints)
Definition: cfModGcd.cc:2030
CanonicalForm monicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition: cfModGcd.cc:2198
long gaussianElimFp(CFMatrix &M, CFArray &L)
Definition: cfModGcd.cc:1739
CanonicalForm cogn
Definition: cfModGcd.cc:4128
int p
Definition: cfModGcd.cc:4077
static CanonicalForm uni_content(const CanonicalForm &F)
compute the content of F, where F is considered as an element of
Definition: cfModGcd.cc:281
CFArray solveGeneralVandermonde(const CFArray &M, const CFArray &A)
Definition: cfModGcd.cc:1632
cl
Definition: cfModGcd.cc:4099
f
Definition: cfModGcd.cc:4080
CanonicalForm extractContents(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &contentF, CanonicalForm &contentG, CanonicalForm &ppF, CanonicalForm &ppG, const int d)
Definition: cfModGcd.cc:311
CanonicalForm lcf
Definition: cfModGcd.cc:4096
static CanonicalForm uni_lcoeff(const CanonicalForm &F)
compute the leading coefficient of F, where F is considered as an element of , order on is dp.
Definition: cfModGcd.cc:339
g
Definition: cfModGcd.cc:4089
int maxNumVars
Definition: cfModGcd.cc:4129
CanonicalForm fp
Definition: cfModGcd.cc:4101
CFArray solveVandermonde(const CFArray &M, const CFArray &A)
Definition: cfModGcd.cc:1579
int i
Definition: cfModGcd.cc:4077
CanonicalForm cogp
Definition: cfModGcd.cc:4128
const CanonicalForm const CanonicalForm & coF
Definition: cfModGcd.cc:67
CanonicalForm cofn
Definition: cfModGcd.cc:4128
CanonicalForm cof
Definition: cfModGcd.cc:4128
const CanonicalForm & GG
Definition: cfModGcd.cc:4075
CanonicalForm cd(bCommonDen(FF))
Definition: cfModGcd.cc:4088
CanonicalForm cog
Definition: cfModGcd.cc:4128
CanonicalForm nonMonicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition: cfModGcd.cc:2482
int minCommonDeg
Definition: cfModGcd.cc:4103
CanonicalForm modGCDFp(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, bool &topLevel, CFList &l)
Definition: cfModGcd.cc:1223
CFArray solveSystemFp(const CFMatrix &M, const CFArray &L)
Definition: cfModGcd.cc:1839
const CanonicalForm & G
Definition: cfModGcd.cc:66
CanonicalForm sparseGCDFp(const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l)
Definition: cfModGcd.cc:3561
static CanonicalForm randomElement(const CanonicalForm &F, const Variable &alpha, CFList &list, bool &fail)
compute a random element a of , s.t. F(a) , F is a univariate polynomial, returns fail if there are...
Definition: cfModGcd.cc:379
CanonicalForm gcdcfcg
Definition: cfModGcd.cc:4100
CanonicalForm cf
Definition: cfModGcd.cc:4082
CanonicalForm Dn
Definition: cfModGcd.cc:4095
CanonicalForm newCof
Definition: cfModGcd.cc:4128
static CanonicalForm FpRandomElement(const CanonicalForm &F, CFList &list, bool &fail)
Definition: cfModGcd.cc:1171
CanonicalForm gp
Definition: cfModGcd.cc:4101
CanonicalForm modGCDGF(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, CFList &l, bool &topLevel)
GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for Computer Algebra" by Gedde...
Definition: cfModGcd.cc:872
static CanonicalForm newtonInterp(const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x)
Newton interpolation - Incremental algorithm. Given a list of values alpha_i and a list of polynomial...
Definition: cfModGcd.cc:364
bool equal
Definition: cfModGcd.cc:4125
CFArray solveSystemFq(const CFMatrix &M, const CFArray &L, const Variable &alpha)
Definition: cfModGcd.cc:1891
CanonicalForm test
Definition: cfModGcd.cc:4095
CanonicalForm b
Definition: cfModGcd.cc:4102
CanonicalForm cg
Definition: cfModGcd.cc:4082
CanonicalForm cDn
Definition: cfModGcd.cc:4128
int d_deg
Definition: cfModGcd.cc:4077
CFList evaluationPoints(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Feval, CanonicalForm &Geval, const CanonicalForm &LCF, const bool &GF, const Variable &alpha, bool &fail, CFList &list)
Definition: cfModGcd.cc:2047
modular and sparse modular GCD algorithms over finite fields and Z.
bool terminationTest(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &coF, const CanonicalForm &coG, const CanonicalForm &cand)
CanonicalForm modGCDZ(const CanonicalForm &FF, const CanonicalForm &GG)
modular GCD over Z
static void evalPoint(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &FEval, CanonicalForm &GEval, CFGenerator &evalPoint)
This file provides functions to compute the Newton polygon of a bivariate polynomial.
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.
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
static const int SW_USE_CHINREM_GCD
set to 1 to use modular gcd over Z
Definition: cf_defs.h:41
#define DELETE_ARRAY(P)
Definition: cf_defs.h:65
#define NEW_ARRAY(T, N)
Definition: cf_defs.h:64
static CanonicalForm balance_p(const CanonicalForm &f, const CanonicalForm &q, const CanonicalForm &qh)
Definition: cf_gcd.cc:149
CanonicalForm randomIrredpoly(int i, const Variable &x)
computes a random monic irreducible univariate polynomial in x over Fp of degree i via NTL/FLINT
Definition: cf_irred.cc:26
generate random irreducible univariate polynomials
Iterators for CanonicalForm's.
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
map polynomials
CanonicalForm mapPrimElem(const CanonicalForm &primElem, const Variable &alpha, const Variable &beta)
compute the image of a primitive element of in . We assume .
Definition: cf_map_ext.cc:450
CanonicalForm GFMapDown(const CanonicalForm &F, int k)
maps a polynomial over to a polynomial over , d needs to be a multiple of k
Definition: cf_map_ext.cc:276
CanonicalForm primitiveElement(const Variable &alpha, Variable &beta, bool &fail)
determine a primitive element of , is a primitive element of a field which is isomorphic to
Definition: cf_map_ext.cc:342
static CanonicalForm mapDown(const CanonicalForm &F, const Variable &alpha, const CanonicalForm &G, CFList &source, CFList &dest)
the CanonicalForm G is the output of map_up, returns F considered as an element over ,...
Definition: cf_map_ext.cc:123
static CanonicalForm mapUp(const Variable &alpha, const Variable &beta)
and is a primitive element, returns the image of
Definition: cf_map_ext.cc:70
CanonicalForm GFMapUp(const CanonicalForm &F, int k)
maps a polynomial over to a polynomial over , d needs to be a multiple of k
Definition: cf_map_ext.cc:240
This file implements functions to map between extensions of finite fields.
int cf_getBigPrime(int i)
Definition: cf_primes.cc:39
int cf_getNumBigPrimes()
Definition: cf_primes.cc:45
access to prime tables
GLOBAL_VAR flint_rand_t FLINTrandom
Definition: cf_random.cc:25
generate random integers, random elements of finite fields
generate random evaluation points
VAR void(* factoryError)(const char *s)
Definition: cf_util.cc:80
int ipower(int b, int m)
int ipower ( int b, int m )
Definition: cf_util.cc:27
generate random elements in F_p(alpha)
Definition: cf_random.h:70
CanonicalForm generate() const
Definition: cf_random.cc:157
int size() const
Definition: ftmpl_array.cc:92
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
CanonicalForm genOne() const
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 inCoeffDomain() const
int level() const
level() returns the level of CO.
bool inBaseDomain() const
CanonicalForm lc() const
CanonicalForm CanonicalForm::lc (), Lc (), LC (), LC ( v ) const.
bool isUnivariate() const
generate random elements in F_p
Definition: cf_random.h:44
CanonicalForm generate() const
Definition: cf_random.cc:68
generate random elements in GF
Definition: cf_random.h:32
CanonicalForm generate() const
Definition: cf_random.cc:78
int length() const
Definition: ftmpl_list.cc:273
void append(const T &)
Definition: ftmpl_list.cc:256
T getLast() const
Definition: ftmpl_list.cc:309
int rows() const
Definition: ftmpl_matrix.h:43
int columns() const
Definition: ftmpl_matrix.h:44
factory's class for variables
Definition: factory.h:127
int level() const
Definition: factory.h:143
functions to print debug output
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
CFFListIterator iter
Definition: facAbsBiFact.cc:53
Variable alpha
Definition: facAbsBiFact.cc:51
return result
Definition: facAbsBiFact.cc:75
CanonicalForm Feval
Definition: facAbsFact.cc:60
CanonicalForm H
Definition: facAbsFact.cc:60
CanonicalForm mipo
Definition: facAlgExt.cc:57
b *CanonicalForm B
Definition: facBivar.cc:52
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
CanonicalForm LCF
Definition: facFactorize.cc:52
CFList & eval
Definition: facFactorize.cc:47
CFList *& Aeval
<[in] poly
Definition: facFactorize.h:31
CFList tmp1
Definition: facFqBivar.cc:72
CFList tmp2
Definition: facFqBivar.cc:72
CFList evalPoints(const CanonicalForm &F, CFList &eval, const Variable &alpha, CFList &list, const bool &GF, bool &fail)
evaluation point search for multivariate factorization, looks for a (F.level() - 1)-tuple such that t...
int j
Definition: facHensel.cc:110
fq_nmod_ctx_clear(fq_con)
nmod_poly_init(FLINTmipo, getCharacteristic())
fq_nmod_ctx_init_modulus(fq_con, FLINTmipo, "Z")
convertFacCF2nmod_poly_t(FLINTmipo, M)
nmod_poly_clear(FLINTmipo)
This file defines functions for Hensel lifting.
Variable FACTORY_PUBLIC rootOf(const CanonicalForm &, char name='@')
returns a symbolic root of polynomial with name name Use it to define algebraic variables
Definition: variable.cc:162
void prune1(const Variable &alpha)
Definition: variable.cc:291
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition: variable.cc:207
void FACTORY_PUBLIC prune(Variable &alpha)
Definition: variable.cc:261
void setMipo(const Variable &alpha, const CanonicalForm &mipo)
Definition: variable.cc:219
#define const
Definition: fegetopt.c:39
some useful template functions.
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)
template bool find(const List< CanonicalForm > &, const CanonicalForm &)
#define D(A)
Definition: gentable.cc:131
VAR char gf_name
Definition: gfops.cc:52
static long ind2(long arg)
Definition: kstd2.cc:532
gmp_float log(const gmp_float &a)
Definition: mpr_complex.cc:343
#define NULL
Definition: omList.c:12
int status int void * buf
Definition: si_signals.h:59
#define A
Definition: sirandom.c:24
#define M
Definition: sirandom.c:25
#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
int gcd(int a, int b)
Definition: walkSupport.cc:836