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