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