source: git/factory/cf_gcd_smallp.cc @ 911444

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