source: git/factory/cf_gcd_smallp.cc @ ad316f8

spielwiese
Last change on this file since ad316f8 was ad316f8, checked in by Hans Schoenemann <hannes@…>, 14 years ago
SunOS port git-svn-id: file:///usr/local/Singular/svn/trunk@12767 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100755
File size: 29.4 KB
Line 
1// -*- c++ -*-
2//*****************************************************************************
3/** @file cf_gcd_smallp.cc
4 *
5 * @author Martin Lee
6 * @date 22.10.2009
7 *
8 * This file implements the GCD of two polynomials over \f$ F_{p} \f$ ,
9 * \f$ F_{p}(\alpha ) \f$ or GF based on Alg. 7.2. as described in "Algorithms
10 * for Computer Algebra" by Geddes, Czapor, Labahnn
11 *
12 * @par Copyright:
13 *   (c) by The SINGULAR Team, see LICENSE file
14 *
15 * @internal
16 * @version \$Id$
17 *
18**/
19//*****************************************************************************
20
21#include <config.h>
22
23#include "assert.h"
24#include "debug.h"
25#include "timing.h"
26
27#include "canonicalform.h"
28#include "cf_map.h"
29#include "ftmpl_functions.h"
30#include "cf_map_ext.cc"
31#include "cf_random.h"
32
33#ifdef HAVE_NTL
34#include <NTL/ZZ_pEX.h>
35#endif
36
37#ifdef HAVE_NTL
38
39TIMING_DEFINE_PRINT(gcd_recursion);
40TIMING_DEFINE_PRINT(newton_interpolation);
41
42/// compressing two polynomials F and G, M is used for compressing,
43/// N to reverse the compression
44int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
45                CFMap & N, bool& top_level) 
46{ 
47  int n= tmax (F.level(), G.level());
48  int * degsf= new int [n + 1];
49  int * degsg= new int [n + 1];
50
51  for (int i = 0; i <= n; i++)
52    degsf[i]= degsg[i]= 0;
53   
54  degsf= degrees (F, degsf);
55  degsg= degrees (G, degsg);
56 
57  int both_non_zero= 0;
58  int f_zero= 0;
59  int g_zero= 0;
60  int both_zero= 0;
61
62  if (top_level) 
63  {
64    for (int i= 1; i <= n; i++) 
65    {
66      if (degsf[i] != 0 && degsg[i] != 0) 
67      {
68        both_non_zero++;
69        continue;
70      }
71      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level()) 
72      {
73        f_zero++;
74        continue;
75      }
76      if (degsg[i] == 0 && degsf[i] && i <= F.level()) 
77      {
78        g_zero++;
79        continue;
80      }
81    }
82
83    if (both_non_zero == 0) return 0;
84
85    // map Variables which do not occur in both polynomials to higher levels
86    int k= 1;
87    int l= 1;
88    for (int i= 1; i <= n; i++) 
89    { 
90      if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level()) 
91      {
92        if (k + both_non_zero != i) 
93        {
94          M.newpair (Variable (i), Variable (k + both_non_zero));
95          N.newpair (Variable (k + both_non_zero), Variable (i));
96        }
97        k++;
98      }
99      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level()) 
100      {
101        if (l + g_zero + both_non_zero != i) 
102        {
103          M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
104          N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
105        }
106        l++;
107      }
108    }
109 
110    // sort Variables x_{i} in increasing order of
111    // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
112    int m= tmin (F.level(), G.level());
113    int max_min_deg;
114    k= both_non_zero;
115    l= 0;
116    int i= 1;
117    while (k > 0) 
118    {
119      max_min_deg= tmin (degsf[i], degsg[i]);
120      while (max_min_deg == 0) 
121      {
122        i++;
123        max_min_deg= tmin (degsf[i], degsg[i]);
124      }
125      for (int j= i + 1; j <=  m; j++) 
126      {
127        if (tmin (degsf[j],degsg[j]) >= max_min_deg) 
128        {
129          max_min_deg= tmin (degsf[j], degsg[j]);
130          l= j; 
131        }
132      }
133      if (l != 0) 
134      {
135        if (l != k) 
136        {
137          M.newpair (Variable (l), Variable(k));
138          N.newpair (Variable (k), Variable(l));
139          degsf[l]= 0;
140          degsg[l]= 0;
141          l= 0;
142        }
143        else 
144        {
145          degsf[l]= 0;
146          degsg[l]= 0;
147          l= 0;
148        }
149      } 
150      else if (l == 0) 
151      {
152        if (i != k) 
153        {
154          M.newpair (Variable (i), Variable (k));
155          N.newpair (Variable (k), Variable (i));
156          degsf[i]= 0;
157          degsg[i]= 0;
158        }
159        else 
160        {
161          degsf[i]= 0;
162          degsg[i]= 0;
163        }
164        i++;
165      } 
166      k--;
167    }
168  }
169  else 
170  {
171    //arrange Variables such that no gaps occur
172    for (int i= 1; i <= n; i++) 
173    {
174      if (degsf[i] == 0 && degsg[i] == 0) 
175      {
176        both_zero++;
177        continue;
178      }
179      else 
180      {
181        if (both_zero != 0) 
182        {
183          M.newpair (Variable (i), Variable (i - both_zero));
184          N.newpair (Variable (i - both_zero), Variable (i));
185        }
186      }
187    }
188  }
189
190  delete [] degsf;
191  delete [] degsg;
192
193  return 1;
194}
195
196/// compute the content of F, where F is considered as an element of
197/// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$
198static inline CanonicalForm
199uni_content (const CanonicalForm & F) 
200{ 
201  if (F.inBaseDomain())
202    return F.genOne();
203  if (F.level() == 1 && F.isUnivariate())
204    return F;
205  if (F.level() != 1 && F.isUnivariate())
206    return F.genOne();
207  if (degree (F,1) == 0) return F.genOne();
208
209  int l= F.level();
210  if (l == 2) 
211    return content(F);
212  else 
213  {
214    CanonicalForm pol, c = 0;
215    CFIterator i = F;
216    for (; i.hasTerms(); i++) 
217    {
218      pol= i.coeff(); 
219      pol= uni_content (pol);
220      c= gcd (c, pol);
221      if (c.isOne())
222        return c;
223    }
224    return c;
225  }
226}
227
228/// compute the leading coefficient of F, where F is considered as an element
229/// of \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$, order on Mon (x_{2},\ldots ,x_{n})
230/// is dp.
231static inline
232CanonicalForm uni_lcoeff (const CanonicalForm& F) 
233{
234  if (F.level() <= 1)
235    return F; 
236  else
237  {
238    Variable x= Variable (2);
239    int deg= totaldegree (F, x, F.mvar());
240    for (CFIterator i= F; i.hasTerms(); i++)
241    {
242      if (i.exp() + totaldegree (i.coeff(), x, i.coeff().mvar()) == deg)
243        return uni_lcoeff (i.coeff()); 
244    }
245  }
246}
247
248/// Newton interpolation - Incremental algorithm.
249/// Given a list of values alpha_i and a list of polynomials u_i, 1 <= i <= n,
250/// computes the interpolation polynomial assuming that
251/// the polynomials in u are the results of evaluating the variabe x
252/// of the unknown polynomial at the alpha values.
253/// This incremental version receives only the values of alpha_n and u_n and
254/// the previous interpolation polynomial for points 1 <= i <= n-1, and computes
255/// the polynomial interpolating in all the points.
256/// newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1})
257static inline CanonicalForm
258newtonInterp(const CanonicalForm alpha, const CanonicalForm u, const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly, const Variable & x)
259{
260  CanonicalForm interPoly;
261
262  interPoly = oldInterPoly + ((u - oldInterPoly(alpha, x)) / newtonPoly(alpha, x)) * newtonPoly;
263  return interPoly;
264}
265
266/// compute a random element a of \f$ F_{p}(\alpha ) \f$ ,
267/// s.t. F(a) \f$ \neq 0 \f$ , F is a univariate polynomial, returns
268/// fail if there are no field elements left which have not been used before
269static inline CanonicalForm
270randomElement (const CanonicalForm & F, const Variable & alpha, CFList & list,
271               bool & fail) 
272{
273  fail= false;
274  Variable x= F.mvar();
275  AlgExtRandomF genAlgExt (alpha);
276  FFRandom genFF;
277  CanonicalForm random, mipo;
278  mipo= getMipo (alpha);
279  int p= getCharacteristic ();
280  int d= degree (mipo);
281  double bound= pow ((double) p, (double) d);
282  do 
283  {
284    if (list.length() == bound)
285    {
286      fail= true;
287      break;
288    }
289    if (list.length() < p) 
290    {
291      random= genFF.generate();
292      while (find (list, random))
293        random= genFF.generate();
294    }
295    else 
296    {
297      random= genAlgExt.generate();
298      while (find (list, random))
299        random= genAlgExt.generate();
300    }
301    if (F (random, x) == 0) 
302    {
303      list.append (random);
304      continue;
305    }
306  } while (find (list, random));
307  return random;
308}
309
310/// chooses a suitable extension of \f$ F_{p}(\alpha ) \f$
311/// we do not extend \f$ F_{p}(\alpha ) \f$ itself but extend \f$ F_{p} \f$ ,
312/// s.t. \f$ F_{p}(\alpha ) \subset F_{p}(\beta ) \f$
313static inline
314void choose_extension (const int& d, const int& num_vars, Variable& beta) 
315{
316  int p= getCharacteristic();
317  ZZ NTLp= to_ZZ (p);
318  ZZ_p::init (NTLp);
319  ZZ_pX NTLirredpoly;
320  //TODO: replace d by max_{i} (deg_x{i}(f))
321  int i= (int) (log ((double) ::pow (d + 1, num_vars))/log ((double) p)); 
322  int m= degree (getMipo (beta));
323  if (i <= 1)
324    i= 2;
325  BuildIrred (NTLirredpoly, i*m); 
326  CanonicalForm mipo= convertNTLZZpX2CF (NTLirredpoly, Variable(1)); 
327  beta= rootOf (mipo); 
328}
329
330/// GCD of F and G over \f$ F_{p}(\alpha ) \f$ ,
331/// l and top_level are only used internally, output is monic
332/// based on Alg. 7.2. as described in "Algorithms for
333/// Computer Algebra" by Geddes, Czapor, Labahn
334CanonicalForm
335GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G, 
336                  Variable & alpha, CFList& l, bool& top_level) 
337{ 
338  CanonicalForm A= F;
339  CanonicalForm B= G;
340  if (F.isZero() && degree(G) > 0) return G/Lc(G);
341  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
342  else if (F.isZero() && G.isZero()) return F.genOne();
343  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
344  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
345  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
346  if (F == G) return F/Lc(F);
347 
348  CFMap M,N;
349  int best_level= myCompress (A, B, M, N, top_level);
350
351  if (best_level == 0) return B.genOne();
352
353  A= M(A);
354  B= M(B);
355
356  Variable x= Variable(1);
357
358  //univariate case
359  if (A.isUnivariate() && B.isUnivariate()) 
360    return N (gcd(A,B)); 
361 
362  CanonicalForm cA, cB;    // content of A and B
363  CanonicalForm ppA, ppB;    // primitive part of A and B
364  CanonicalForm gcdcAcB;
365
366  cA = uni_content (A);
367  cB = uni_content (B);
368
369  gcdcAcB= gcd (cA, cB);
370
371  ppA= A/cA;
372  ppB= B/cB;
373
374  CanonicalForm lcA, lcB;  // leading coefficients of A and B
375  CanonicalForm gcdlcAlcB; 
376
377  lcA= uni_lcoeff (ppA);
378  lcB= uni_lcoeff (ppB);
379 
380  if (fdivides (lcA, lcB)) { 
381    if (fdivides (A, B))
382      return F/Lc(F);   
383  }
384  if (fdivides (lcB, lcA))
385  { 
386    if (fdivides (B, A)) 
387      return G/Lc(G);
388  }
389
390  gcdlcAlcB= gcd (lcA, lcB);
391
392  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
393
394  if (d == 0) return N(gcdcAcB);
395  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
396  if (d0 < d) 
397    d= d0; 
398  if (d == 0) return N(gcdcAcB);
399
400  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
401  CanonicalForm newtonPoly;
402
403  newtonPoly= 1;
404  m= gcdlcAlcB;
405  G_m= 0;
406  H= 0;
407  bool fail= false;
408  top_level= false;
409  bool inextension= false;
410  Variable V_buf= alpha;
411  CanonicalForm prim_elem, im_prim_elem;
412  CFList source, dest;
413  do 
414  {
415    random_element= randomElement (m, V_buf, l, fail);
416    if (fail) 
417    {
418      source= CFList();
419      dest= CFList();
420      int num_vars= tmin (getNumVars(A), getNumVars(B));
421      num_vars--;
422
423      choose_extension (d, num_vars, V_buf);
424      bool prim_fail= false;
425      Variable V_buf2;
426      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
427
428      ASSERT (!prim_fail, "failure in integer factorizer");
429      if (prim_fail)
430        ; //ERROR
431      else
432        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
433
434      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
435      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
436      inextension= true;
437      for (CFListIterator i= l; i.hasItem(); i++) 
438        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem, 
439                             im_prim_elem, source, dest);
440      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
441      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
442      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, 
443                          source, dest);
444      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
445      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
446      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, 
447                         source, dest);
448
449      fail= false;
450      random_element= randomElement (m, V_buf, l, fail );
451      DEBOUTLN (cerr, "fail= " << fail);
452      CFList list;
453      TIMING_START (gcd_recursion);
454      G_random_element= 
455      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
456                        list, top_level);
457      TIMING_END_AND_PRINT (gcd_recursion, 
458                            "time for recursive call: ");
459      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
460    }
461    else 
462    {
463      CFList list;
464      TIMING_START (gcd_recursion);
465      G_random_element= 
466      GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf,
467                        list, top_level);
468      TIMING_END_AND_PRINT (gcd_recursion, 
469                            "time for recursive call: ");
470      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
471    }
472
473    d0= totaldegree (G_random_element, Variable(2), 
474                     Variable (G_random_element.level()));
475    if (d0 == 0) return N(gcdcAcB); 
476    if (d0 >  d) 
477    {
478      if (!find (l, random_element))
479        l.append (random_element);
480      continue;
481    }
482
483    G_random_element= 
484    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
485    * G_random_element;
486
487    d0= totaldegree (G_random_element, Variable(2), 
488                     Variable(G_random_element.level()));
489    if (d0 <  d) 
490    {
491      m= gcdlcAlcB;
492      newtonPoly= 1;
493      G_m= 0;
494      d= d0;
495    }
496
497    TIMING_START (newton_interpolation);
498    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
499    TIMING_END_AND_PRINT (newton_interpolation, 
500                          "time for newton interpolation: ");
501
502    //termination test
503    if (uni_lcoeff (H) == gcdlcAlcB) 
504    {
505      cH= uni_content (H);
506      ppH= H/cH;
507      if (inextension) 
508      {
509        CFList u, v; 
510        //maybe it's better to test if ppH is an element of F(\alpha) before
511        //mapping down
512        DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
513        ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v); 
514        ppH /= Lc(ppH);
515        DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
516        if (fdivides (ppH, A) && fdivides (ppH, B))         
517          return N(gcdcAcB*ppH);
518      }
519      else if (fdivides (ppH, A) && fdivides (ppH, B)) 
520        return N(gcdcAcB*(ppH/Lc(ppH)));
521    }
522
523    G_m= H;
524    newtonPoly= newtonPoly*(x - random_element);
525    m= m*(x - random_element);
526    if (!find (l, random_element))
527      l.append (random_element);
528  } while (1);
529}
530
531/// compute a random element a of GF, s.t. F(a) \f$ \neq 0 \f$ , F is a
532/// univariate polynomial, returns fail if there are no field elements left
533/// which have not been used before
534CanonicalForm
535GFRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
536{
537  fail= false;
538  Variable x= F.mvar();
539  GFRandom genGF;
540  CanonicalForm random;
541  int p= getCharacteristic();
542  int d= getGFDegree();
543  double bound= pow ((double) p, (double) d);
544  do 
545  {
546    if (list.length() == bound)
547    {
548      fail= true;
549      break;
550    }
551    if (list.length() < 1)
552      random= 0;
553    else 
554    {
555      random= genGF.generate();
556      while (find (list, random))
557        random= genGF.generate();
558    }
559    if (F (random, x) == 0) 
560    {
561      list.append (random);
562      continue;
563    }
564  } while (find (list, random));
565  return random;
566}
567
568/// GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for
569/// Computer Algebra" by Geddes, Czapor, Labahn
570/// Usually this algorithm will be faster than GCD_Fp_extension since GF has
571/// faster field arithmetics, however it might fail if the input is large since
572/// the size of the base field is bounded by 2^16, output is monic
573CanonicalForm
574GCD_GF (const CanonicalForm& F, const CanonicalForm& G, CFList& l, 
575        bool& top_level) 
576{ 
577  CanonicalForm A= F;
578  CanonicalForm B= G;
579  if (F.isZero() && degree(G) > 0) return G/Lc(G);
580  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
581  else if (F.isZero() && G.isZero()) return F.genOne();
582  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
583  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
584  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
585  if (F == G) return F/Lc(F);
586 
587  CFMap M,N;
588  int best_level= myCompress (A, B, M, N, top_level);
589
590  if (best_level == 0) return B.genOne();
591
592  A= M(A);
593  B= M(B);
594
595  Variable x= Variable(1);
596
597  //univariate case
598  if (A.isUnivariate() && B.isUnivariate()) 
599    return N (gcd (A, B)); 
600
601  CanonicalForm cA, cB;    // content of A and B
602  CanonicalForm ppA, ppB;    // primitive part of A and B
603  CanonicalForm gcdcAcB;
604
605  cA = uni_content (A);
606  cB = uni_content (B);
607
608  gcdcAcB= gcd (cA, cB);
609
610  ppA= A/cA;
611  ppB= B/cB;
612
613  CanonicalForm lcA, lcB;  // leading coefficients of A and B
614  CanonicalForm gcdlcAlcB; 
615
616  lcA= uni_lcoeff (ppA);
617  lcB= uni_lcoeff (ppB);
618
619  if (fdivides (lcA, lcB)) 
620  { 
621    if (fdivides (A, B))
622      return F; 
623  }   
624  if (fdivides (lcB, lcA)) 
625  { 
626    if (fdivides (B, A)) 
627      return G;
628  }
629
630  gcdlcAlcB= gcd (lcA, lcB);
631
632  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
633  if (d == 0) return N(gcdcAcB);
634  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
635  if (d0 < d) 
636    d= d0; 
637  if (d == 0) return N(gcdcAcB);
638
639  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
640  CanonicalForm newtonPoly;
641
642  newtonPoly= 1;
643  m= gcdlcAlcB;
644  G_m= 0;
645  H= 0;
646  bool fail= false;
647  top_level= false;
648  bool inextension= false;
649  int p;
650  int k= getGFDegree();
651  int kk;
652  double expon;
653  char gf_name_buf= gf_name;
654  do 
655  {
656    random_element= GFRandomElement (m, l, fail);
657    if (fail) 
658    { 
659      int num_vars= tmin (getNumVars(A), getNumVars(B));
660      num_vars--;
661      p= getCharacteristic();
662      expon= floor ((log ((double) ::pow (d + 1, num_vars))/log ((double) p)));
663      if (expon < 2)
664        expon= 2;
665      kk= getGFDegree(); 
666      if (::pow (p, kk*expon) < (1 << 16)) 
667        setCharacteristic (p, kk*(int)expon, 'b');
668      else 
669      {
670        expon= floor((log ((double)((1<<16) - 1)))/(log((double)p)*kk));
671        ASSERT (expon >= 2, "not enough points in GCD_GF");
672        setCharacteristic (p, (int)(kk*expon), 'b');
673      }
674      inextension= true;
675      fail= false;
676      for (CFListIterator i= l; i.hasItem(); i++) 
677        i.getItem()= GFMapUp (i.getItem(), kk);
678      m= GFMapUp (m, kk);
679      G_m= GFMapUp (G_m, kk);
680      newtonPoly= GFMapUp (newtonPoly, kk);
681      ppA= GFMapUp (ppA, kk);
682      ppB= GFMapUp (ppB, kk);
683      gcdlcAlcB= GFMapUp (gcdlcAlcB, kk);
684      random_element= GFRandomElement (m, l, fail);
685      DEBOUTLN (cerr, "fail= " << fail);
686      CFList list;
687      TIMING_START (gcd_recursion);
688      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 
689                                list, top_level);
690      TIMING_END_AND_PRINT (gcd_recursion, 
691                            "time for recursive call: ");
692      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
693    }
694    else 
695    {
696      CFList list;
697      TIMING_START (gcd_recursion);
698      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 
699                                list, top_level);
700      TIMING_END_AND_PRINT (gcd_recursion, 
701                            "time for recursive call: ");
702      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
703    }
704
705    d0= totaldegree (G_random_element, Variable(2), 
706                     Variable (G_random_element.level()));
707    if (d0 == 0) 
708    {
709      if (inextension) 
710      {
711        setCharacteristic (p, k, gf_name_buf);
712        return N(gcdcAcB); 
713      }
714      else
715        return N(gcdcAcB);
716    } 
717    if (d0 >  d) 
718    {
719      if (!find (l, random_element))
720        l.append (random_element);
721      continue;
722    }
723
724    G_random_element= 
725    (gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) *
726     G_random_element;
727    d0= totaldegree (G_random_element, Variable(2), 
728                     Variable (G_random_element.level()));
729
730    if (d0 < d) 
731    {
732      m= gcdlcAlcB;
733      newtonPoly= 1;
734      G_m= 0;
735      d= d0;
736    }
737
738    TIMING_START (newton_interpolation);
739    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
740    TIMING_END_AND_PRINT (newton_interpolation, "time for newton interpolation: ");
741
742    //termination test
743    if (uni_lcoeff (H) == gcdlcAlcB) 
744    {
745      cH= uni_content (H);
746      ppH= H/cH;
747      if (inextension) 
748      {
749        if (fdivides(ppH, GFMapUp(A, k)) && fdivides(ppH, GFMapUp(B,k))) 
750        {
751          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
752          ppH= GFMapDown (ppH, k);
753          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
754          setCharacteristic (p, k, gf_name_buf);
755          return N(gcdcAcB*ppH);
756        }
757      }
758      else 
759      {
760        if (fdivides (ppH, A) && fdivides (ppH, B)) 
761          return N(gcdcAcB*ppH);
762      }
763    }
764
765    G_m= H;
766    newtonPoly= newtonPoly*(x - random_element);
767    m= m*(x - random_element);
768    if (!find (l, random_element))
769      l.append (random_element);
770  } while (1);
771}
772
773/// F is assumed to be an univariate polynomial in x,
774/// computes a random monic irreducible univariate polynomial of random
775/// degree < i in x which does not divide F
776CanonicalForm
777randomIrredpoly (int i, const Variable & x, 
778                 CFList & list) 
779{
780  int random;
781  int p= getCharacteristic();
782  ZZ NTLp= to_ZZ (p);
783  ZZ_p::init (NTLp);
784  ZZ_pX NTLirredpoly; 
785  CanonicalForm CFirredpoly;
786  double buf= ceil((log((double) i) / log((double) p)));
787  buf= buf + 1;
788  double bound= 0;
789  //lower bound on the number of monic irreducibles of degree less than buf
790  for (int j= 2; j < buf; j++) 
791    bound += ((pow ((double)p, (double) j) - 2*pow((double) p, 
792              (double) (j/2))) / j);
793  if (list.length() > bound) 
794  {
795    if (buf > 1)
796      buf--;
797    buf *= 2;
798    buf++;
799  }
800  for (int j= ((int) (buf - 1)/2) + 1; j < buf; j++) 
801    bound += ((pow ((double)p, (double) j) - 2*pow((double) p, 
802              (double) (j/2))) / j);
803  do 
804  {
805    if (list.length() <= bound) 
806    {
807      random= factoryrandom ((int) buf); 
808      if (random <= 1)
809        random= 2;
810    }
811    else 
812    {
813      if (buf > 1) 
814        buf--;
815      buf *= 2;
816      buf++;
817      for (int j= ((int) (buf - 1)/2) + 1; j < buf; j++) 
818        bound += ((pow ((double)p, (double) j) - 2*pow((double) p, 
819                  (double) (j/2))) / j);
820      random= factoryrandom ((int) buf);
821      if (random <= 1)
822        random= 2;
823    } 
824    BuildIrred (NTLirredpoly, random);
825    CFirredpoly= convertNTLZZpX2CF (NTLirredpoly, x);
826  } while (find (list, CFirredpoly));
827  return CFirredpoly;
828} 
829
830CanonicalForm
831FpRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
832{
833  fail= false;
834  Variable x= F.mvar();
835  FFRandom genFF;
836  CanonicalForm random;
837  int p= getCharacteristic();
838  double bound= p;
839  do 
840  {
841    if (list.length() == bound)
842    {
843      fail= true;
844      break;
845    }
846    if (list.length() < 1)
847      random= 0;
848    else 
849    {
850      random= genFF.generate();
851      while (find (list, random))
852        random= genFF.generate();
853    }
854    if (F (random, x) == 0) 
855    {
856      list.append (random);
857      continue;
858    }
859  } while (find (list, random));
860  return random;
861}
862
863CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G, 
864                           bool& top_level, CFList& l) 
865{
866  CanonicalForm A= F;
867  CanonicalForm B= G;
868  if (F.isZero() && degree(G) > 0) return G/Lc(G);
869  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
870  else if (F.isZero() && G.isZero()) return F.genOne();
871  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
872  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
873  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
874  if (F == G) return F/Lc(F);
875
876  CFMap M,N;
877  int best_level= myCompress (A, B, M, N, top_level);
878
879  if (best_level == 0) return B.genOne();
880
881  A= M(A);
882  B= M(B);
883
884  Variable x= Variable (1);
885
886  //univariate case
887  if (A.isUnivariate() && B.isUnivariate()) 
888    return N (gcd (A, B));
889
890  CanonicalForm cA, cB;    // content of A and B
891  CanonicalForm ppA, ppB;    // primitive part of A and B
892  CanonicalForm gcdcAcB;
893  cA = uni_content (A);
894  cB = uni_content (B); 
895  gcdcAcB= gcd (cA, cB);
896  ppA= A/cA;
897  ppB= B/cB;
898
899  CanonicalForm lcA, lcB;  // leading coefficients of A and B
900  CanonicalForm gcdlcAlcB; 
901  lcA= uni_lcoeff (ppA);
902  lcB= uni_lcoeff (ppB);
903
904  if (fdivides (lcA, lcB)) 
905  { 
906    if (fdivides (A, B))
907      return F/Lc(F);
908  }   
909  if (fdivides (lcB, lcA))
910  { 
911    if (fdivides (B, A)) 
912      return G/Lc(G);
913  }
914
915  gcdlcAlcB= gcd (lcA, lcB); 
916 
917  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
918  int d0;
919
920  if (d == 0) return N(gcdcAcB);
921  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
922
923  if (d0 < d) 
924    d= d0;
925
926  if (d == 0) return N(gcdcAcB);
927
928  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
929  CanonicalForm newtonPoly= 1;
930  m= gcdlcAlcB;
931  H= 0;
932  G_m= 0;
933  Variable alpha, V_buf;
934  double expon;
935  int p, k;
936  bool fail= false;
937  bool inextension= false;
938  bool inextensionextension= false;
939  top_level= false;
940  CanonicalForm CF_buf;
941  CFList source, dest;
942  CanonicalForm gcdcheck;
943  do 
944  {
945    if (inextension)
946      random_element= randomElement (m, alpha, l, fail);
947    else
948      random_element= FpRandomElement (m, l, fail);
949
950    if (!fail && !inextension)
951    {
952      CFList list;
953      TIMING_START (gcd_recursion);
954      G_random_element= 
955      GCD_small_p (ppA (random_element,x), ppB (random_element,x), top_level, 
956      list);
957      TIMING_END_AND_PRINT (gcd_recursion, 
958                            "time for recursive call: ");
959      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
960    }
961    else if (!fail && inextension)
962    {
963      CFList list; 
964      TIMING_START (gcd_recursion);
965      G_random_element= 
966      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 
967                        list, top_level);
968      TIMING_END_AND_PRINT (gcd_recursion, 
969                            "time for recursive call: ");
970      DEBOUTLN (cerr, "G_random_element= " << G_random_element);   
971    }
972    else if (fail && !inextension)
973    {
974      source= CFList();
975      dest= CFList();
976      CFList list;
977      CanonicalForm mipo;
978      int deg= 3;
979      do {
980        mipo= randomIrredpoly (deg, x, list); 
981        alpha= rootOf (mipo);
982        inextension= true;
983        fail= false;
984        random_element= randomElement (m, alpha, l, fail); 
985        deg++;
986      } while (fail); 
987      list= CFList();
988      TIMING_START (gcd_recursion);
989      G_random_element= 
990      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 
991                        list, top_level);
992      TIMING_END_AND_PRINT (gcd_recursion, 
993                            "time for recursive call: ");
994      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
995    }
996    else if (fail && inextension)
997    {
998      source= CFList();
999      dest= CFList();
1000      int num_vars= tmin (getNumVars(A), getNumVars(B));
1001      num_vars--;
1002      V_buf= alpha;
1003      choose_extension (d, num_vars, V_buf);
1004      bool prim_fail= false;
1005      Variable V_buf2;
1006      CanonicalForm prim_elem, im_prim_elem;
1007      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
1008       
1009      ASSERT (!prim_fail, "failure in integer factorizer");
1010      if (prim_fail)
1011        ; //ERROR
1012      else
1013        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
1014
1015      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
1016      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
1017
1018      inextensionextension= true;
1019      for (CFListIterator i= l; i.hasItem(); i++) 
1020        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem, 
1021                             im_prim_elem, source, dest);
1022      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1023      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1024      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, 
1025                          source, dest);
1026      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1027      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1028      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, 
1029                         source, dest);
1030      fail= false;
1031      random_element= randomElement (m, V_buf, l, fail );
1032      DEBOUTLN (cerr, "fail= " << fail);
1033      CFList list;
1034      TIMING_START (gcd_recursion);
1035      G_random_element= 
1036      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
1037                        list, top_level);
1038      TIMING_END_AND_PRINT (gcd_recursion, 
1039                            "time for recursive call: ");
1040      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1041      alpha= V_buf;
1042    } 
1043
1044    d0= totaldegree (G_random_element, Variable(2), 
1045                     Variable (G_random_element.level()));
1046
1047    if (d0 == 0)
1048    {
1049      return N(gcdcAcB); 
1050    }
1051    if (d0 >  d) 
1052    { 
1053      if (!find (l, random_element))
1054        l.append (random_element);
1055      continue;
1056    }
1057
1058    G_random_element= (gcdlcAlcB(random_element,x)/uni_lcoeff(G_random_element))
1059                       *G_random_element; 
1060
1061   
1062    d0= totaldegree (G_random_element, Variable(2), 
1063                     Variable(G_random_element.level()));
1064
1065    if (d0 <  d) 
1066    {
1067      m= gcdlcAlcB;
1068      newtonPoly= 1;
1069      G_m= 0;
1070      d= d0;
1071    }
1072   
1073    TIMING_START (newton_interpolation);
1074    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1075    TIMING_END_AND_PRINT (newton_interpolation, 
1076                          "time for newton_interpolation: ");
1077
1078    //termination test
1079    if (uni_lcoeff (H) == gcdlcAlcB) 
1080    {
1081      cH= uni_content (H);
1082      ppH= H/cH;
1083      ppH /= Lc (ppH);
1084      DEBOUTLN (cerr, "ppH= " << ppH);
1085      if (fdivides (ppH, A) && fdivides (ppH, B)) 
1086        return N(gcdcAcB*ppH);
1087    }
1088
1089    G_m= H;
1090    newtonPoly= newtonPoly*(x - random_element);
1091    m= m*(x - random_element);
1092    if (!find (l, random_element))
1093      l.append (random_element);
1094  } while (1);
1095}
1096
1097#endif
Note: See TracBrowser for help on using the repository browser.