source: git/factory/cf_gcd_smallp.cc @ 51615d6

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