source: git/factory/cf_gcd_smallp.cc @ 2f573cc

fieker-DuValspielwiese
Last change on this file since 2f573cc was dec1024, checked in by Martin Lee <martinlee84@…>, 14 years ago
added heuristics to modular gcd git-svn-id: file:///usr/local/Singular/svn/trunk@13230 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100755
File size: 34.8 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& topLevel) 
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 (topLevel) 
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) 
91    {
92      delete [] degsf;
93      delete [] degsg;
94      return 0;
95    }
96
97    // map Variables which do not occur in both polynomials to higher levels
98    int k= 1;
99    int l= 1;
100    for (int i= 1; i <= n; i++) 
101    { 
102      if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level()) 
103      {
104        if (k + both_non_zero != i) 
105        {
106          M.newpair (Variable (i), Variable (k + both_non_zero));
107          N.newpair (Variable (k + both_non_zero), Variable (i));
108        }
109        k++;
110      }
111      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level()) 
112      {
113        if (l + g_zero + both_non_zero != i) 
114        {
115          M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
116          N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
117        }
118        l++;
119      }
120    }
121 
122    // sort Variables x_{i} in increasing order of
123    // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
124    int m= tmin (F.level(), G.level());
125    int max_min_deg;
126    k= both_non_zero;
127    l= 0;
128    int i= 1;
129    while (k > 0) 
130    {
131      max_min_deg= tmin (degsf[i], degsg[i]);
132      while (max_min_deg == 0) 
133      {
134        i++;
135        max_min_deg= tmin (degsf[i], degsg[i]);
136      }
137      for (int j= i + 1; j <=  m; j++) 
138      {
139        if (tmin (degsf[j],degsg[j]) >= max_min_deg) 
140        {
141          max_min_deg= tmin (degsf[j], degsg[j]);
142          l= j; 
143        }
144      }
145      if (l != 0) 
146      {
147        if (l != k) 
148        {
149          M.newpair (Variable (l), Variable(k));
150          N.newpair (Variable (k), Variable(l));
151          degsf[l]= 0;
152          degsg[l]= 0;
153          l= 0;
154        }
155        else 
156        {
157          degsf[l]= 0;
158          degsg[l]= 0;
159          l= 0;
160        }
161      } 
162      else if (l == 0) 
163      {
164        if (i != k) 
165        {
166          M.newpair (Variable (i), Variable (k));
167          N.newpair (Variable (k), Variable (i));
168          degsf[i]= 0;
169          degsg[i]= 0;
170        }
171        else 
172        {
173          degsf[i]= 0;
174          degsg[i]= 0;
175        }
176        i++;
177      } 
178      k--;
179    }
180  }
181  else 
182  {
183    //arrange Variables such that no gaps occur
184    for (int i= 1; i <= n; i++) 
185    {
186      if (degsf[i] == 0 && degsg[i] == 0) 
187      {
188        both_zero++;
189        continue;
190      }
191      else 
192      {
193        if (both_zero != 0) 
194        {
195          M.newpair (Variable (i), Variable (i - both_zero));
196          N.newpair (Variable (i - both_zero), Variable (i));
197        }
198      }
199    }
200  }
201
202  delete [] degsf;
203  delete [] degsg;
204
205  return 1;
206}
207
208int
209substituteCheck (const CanonicalForm& F, const CanonicalForm& G) 
210{
211  if (F.inCoeffDomain() || G.inCoeffDomain())
212    return 0;
213  Variable x= Variable (1);
214  if (degree (F, x) <= 1 || degree (G, x) <= 1)
215    return 0;
216  CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping seems to be pretty expensive
217  CanonicalForm g= swapvar (G, G.mvar(), x); 
218  int sizef= 0;
219  int sizeg= 0; 
220  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
221  {
222    if (i.exp() == 1)
223      return 0;
224  }
225  for (CFIterator i= g; i.hasTerms(); i++, sizeg++)
226  {
227    if (i.exp() == 1)
228      return 0;
229  }
230  int * expf= new int [sizef];
231  int * expg= new int [sizeg];
232  int j= 0;
233  for (CFIterator i= f; i.hasTerms(); i++, j++)
234  {
235    expf [j]= i.exp();
236  }
237  j= 0;
238  for (CFIterator i= g; i.hasTerms(); i++, j++)
239  {
240    expg [j]= i.exp();
241  }
242 
243  int indf= sizef - 1;
244  int indg= sizeg - 1;
245  if (expf[indf] == 0)
246    indf--;
247  if (expg[indg] == 0)
248    indg--;
249   
250  if (expg[indg] != expf [indf] || (expg[indg] == 1 && expf[indf] == 1))
251  {
252    delete [] expg;
253    delete [] expf;
254    return 0;
255  }
256  // expf[indg] == expf[indf] after here
257  int result= expg[indg];
258  for (int i= indf - 1; i >= 0; i--)
259  {
260    if (expf [i]%result != 0)
261    {
262      delete [] expg;
263      delete [] expf;
264      return 0;
265    }
266  }
267 
268  for (int i= indg - 1; i >= 0; i--)
269  {
270    if (expg [i]%result != 0)
271    {
272      delete [] expg;
273      delete [] expf;
274      return 0;
275    }
276  }
277
278  delete [] expg;
279  delete [] expf;
280  return result;
281}
282
283// substiution
284void 
285subst (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& A, 
286       CanonicalForm& B, const int d
287      )
288{
289  if (d == 1)
290  {
291    A= F;
292    B= G;
293    return;
294  }
295 
296  CanonicalForm C= 0;
297  CanonicalForm D= 0; 
298  Variable x= Variable (1);
299  CanonicalForm f= swapvar (F, x, F.mvar());
300  CanonicalForm g= swapvar (G, x, G.mvar());
301  for (CFIterator i= f; i.hasTerms(); i++)
302    C += i.coeff()*power (f.mvar(), i.exp()/ d);
303  for (CFIterator i= g; i.hasTerms(); i++)
304    D += i.coeff()*power (g.mvar(), i.exp()/ d);
305  A= swapvar (C, x, F.mvar());
306  B= swapvar (D, x, G.mvar());
307}
308
309CanonicalForm
310reverseSubst (const CanonicalForm& F, const int d)
311{
312  if (d == 1)
313    return F;
314 
315  Variable x= Variable (1);
316  if (degree (F, x) <= 0)
317    return F;
318  CanonicalForm f= swapvar (F, x, F.mvar());
319  CanonicalForm result= 0;
320  for (CFIterator i= f; i.hasTerms(); i++)
321    result += i.coeff()*power (f.mvar(), d*i.exp());
322  return swapvar (result, x, F.mvar());   
323}
324
325static inline CanonicalForm
326uni_content (const CanonicalForm & F);
327
328CanonicalForm
329uni_content (const CanonicalForm& F, const Variable& x)
330{
331  if (F.inCoeffDomain())
332    return F.genOne();
333  if (F.level() == x.level() && F.isUnivariate())
334    return F;
335  if (F.level() != x.level() && F.isUnivariate())
336    return F.genOne();
337   
338  if (x.level() != 1)
339  {
340    CanonicalForm f= swapvar (F, x, Variable (1));
341    CanonicalForm result= uni_content (f);
342    return swapvar (result, x, Variable (1));
343  }
344  else
345    return uni_content (F);
346} 
347
348/// compute the content of F, where F is considered as an element of
349/// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$
350static inline CanonicalForm
351uni_content (const CanonicalForm & F) 
352{ 
353  if (F.inBaseDomain())
354    return F.genOne();
355  if (F.level() == 1 && F.isUnivariate())
356    return F;
357  if (F.level() != 1 && F.isUnivariate())
358    return F.genOne();
359  if (degree (F,1) == 0) return F.genOne();
360
361  int l= F.level();
362  if (l == 2) 
363    return content(F);
364  else 
365  {
366    CanonicalForm pol, c = 0;
367    CFIterator i = F;
368    for (; i.hasTerms(); i++) 
369    {
370      pol= i.coeff(); 
371      pol= uni_content (pol);
372      c= gcd (c, pol);
373      if (c.isOne())
374        return c;
375    }
376    return c;
377  }
378}
379
380CanonicalForm
381extractContents (const CanonicalForm& F, const CanonicalForm& G, 
382                 CanonicalForm& contentF, CanonicalForm& contentG, 
383                 CanonicalForm& ppF, CanonicalForm& ppG, const int d)
384{
385  CanonicalForm uniContentF, uniContentG, gcdcFcG;
386  contentF= 1;
387  contentG= 1;
388  ppF= F;
389  ppG= G;
390  CanonicalForm result= 1;
391  for (int i= 1; i <= d; i++)
392  {
393    uniContentF= uni_content (F, Variable (i));
394    uniContentG= uni_content (G, Variable (i));
395    gcdcFcG= gcd (uniContentF, uniContentG);
396    contentF *= uniContentF;
397    contentG *= uniContentG;
398    ppF /= uniContentF;
399    ppG /= uniContentG;
400    result *= gcdcFcG;
401  }
402  return result;
403}
404
405/// compute the leading coefficient of F, where F is considered as an element
406/// of \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$, order on Mon (x_{2},\ldots ,x_{n})
407/// is dp.
408static inline
409CanonicalForm uni_lcoeff (const CanonicalForm& F) 
410{
411  if (F.level() <= 1)
412    return F; 
413  else
414  {
415    Variable x= Variable (2);
416    int deg= totaldegree (F, x, F.mvar());
417    for (CFIterator i= F; i.hasTerms(); i++)
418    {
419      if (i.exp() + totaldegree (i.coeff(), x, i.coeff().mvar()) == deg)
420        return uni_lcoeff (i.coeff()); 
421    }
422  }
423}
424
425/// Newton interpolation - Incremental algorithm.
426/// Given a list of values alpha_i and a list of polynomials u_i, 1 <= i <= n,
427/// computes the interpolation polynomial assuming that
428/// the polynomials in u are the results of evaluating the variabe x
429/// of the unknown polynomial at the alpha values.
430/// This incremental version receives only the values of alpha_n and u_n and
431/// the previous interpolation polynomial for points 1 <= i <= n-1, and computes
432/// the polynomial interpolating in all the points.
433/// newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1})
434static inline CanonicalForm
435newtonInterp(const CanonicalForm alpha, const CanonicalForm u, const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly, const Variable & x)
436{
437  CanonicalForm interPoly;
438
439  interPoly = oldInterPoly + ((u - oldInterPoly(alpha, x)) / newtonPoly(alpha, x)) * newtonPoly;
440  return interPoly;
441}
442
443/// compute a random element a of \f$ F_{p}(\alpha ) \f$ ,
444/// s.t. F(a) \f$ \neq 0 \f$ , F is a univariate polynomial, returns
445/// fail if there are no field elements left which have not been used before
446static inline CanonicalForm
447randomElement (const CanonicalForm & F, const Variable & alpha, CFList & list,
448               bool & fail) 
449{
450  fail= false;
451  Variable x= F.mvar();
452  AlgExtRandomF genAlgExt (alpha);
453  FFRandom genFF;
454  CanonicalForm random, mipo;
455  mipo= getMipo (alpha);
456  int p= getCharacteristic ();
457  int d= degree (mipo);
458  int bound= ipower (p, d);
459  do 
460  {
461    if (list.length() == bound)
462    {
463      fail= true;
464      break;
465    }
466    if (list.length() < p) 
467    {
468      random= genFF.generate();
469      while (find (list, random))
470        random= genFF.generate();
471    }
472    else 
473    {
474      random= genAlgExt.generate();
475      while (find (list, random))
476        random= genAlgExt.generate();
477    }
478    if (F (random, x) == 0) 
479    {
480      list.append (random);
481      continue;
482    }
483  } while (find (list, random));
484  return random;
485}
486
487/// chooses a suitable extension of \f$ F_{p}(\alpha ) \f$
488/// we do not extend \f$ F_{p}(\alpha ) \f$ itself but extend \f$ F_{p} \f$ ,
489/// s.t. \f$ F_{p}(\alpha ) \subset F_{p}(\beta ) \f$
490static inline
491void choose_extension (const int& d, const int& num_vars, Variable& beta) 
492{
493  int p= getCharacteristic();
494  ZZ NTLp= to_ZZ (p);
495  ZZ_p::init (NTLp);
496  ZZ_pX NTLirredpoly;
497  //TODO: replace d by max_{i} (deg_x{i}(f))
498  int i= (int) (log ((double) ipower (d + 1, num_vars))/log ((double) p)); 
499  int m= degree (getMipo (beta));
500  if (i <= 1)
501    i= 2;
502  BuildIrred (NTLirredpoly, i*m); 
503  CanonicalForm mipo= convertNTLZZpX2CF (NTLirredpoly, Variable(1)); 
504  beta= rootOf (mipo); 
505}
506
507/// GCD of F and G over \f$ F_{p}(\alpha ) \f$ ,
508/// l and topLevel are only used internally, output is monic
509/// based on Alg. 7.2. as described in "Algorithms for
510/// Computer Algebra" by Geddes, Czapor, Labahn
511CanonicalForm
512GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G, 
513                  Variable & alpha, CFList& l, bool& topLevel) 
514{ 
515  CanonicalForm A= F;
516  CanonicalForm B= G;
517  if (F.isZero() && degree(G) > 0) return G/Lc(G);
518  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
519  else if (F.isZero() && G.isZero()) return F.genOne();
520  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
521  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
522  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
523  if (F == G) return F/Lc(F);
524 
525  CFMap M,N;
526  int best_level= myCompress (A, B, M, N, topLevel);
527
528  if (best_level == 0) return B.genOne();
529
530  A= M(A);
531  B= M(B);
532
533  Variable x= Variable(1);
534
535  //univariate case
536  if (A.isUnivariate() && B.isUnivariate()) 
537    return N (gcd(A,B)); 
538 
539  int substitute= substituteCheck (A, B);
540 
541  if (substitute > 1)
542    subst (A, B, A, B, substitute);
543
544  CanonicalForm cA, cB;    // content of A and B
545  CanonicalForm ppA, ppB;    // primitive part of A and B
546  CanonicalForm gcdcAcB;
547
548  if (topLevel)
549  {
550    if (best_level <= 2)
551      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
552    else 
553      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 
554  }
555  else
556  {
557    cA = uni_content (A);
558    cB = uni_content (B); 
559    gcdcAcB= gcd (cA, cB);
560    ppA= A/cA;
561    ppB= B/cB;
562  }
563
564  CanonicalForm lcA, lcB;  // leading coefficients of A and B
565  CanonicalForm gcdlcAlcB; 
566
567  lcA= uni_lcoeff (ppA);
568  lcB= uni_lcoeff (ppB);
569 
570  if (fdivides (lcA, lcB)) 
571  { 
572    if (fdivides (A, B))
573      return F/Lc(F);   
574  }
575  if (fdivides (lcB, lcA))
576  { 
577    if (fdivides (B, A)) 
578      return G/Lc(G);
579  }
580
581  gcdlcAlcB= gcd (lcA, lcB);
582
583  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
584
585  if (d == 0)
586  {
587    if (substitute > 1)
588      return N (reverseSubst (gcdcAcB, substitute));
589    else 
590      return N(gcdcAcB);
591  }
592  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
593  if (d0 < d) 
594    d= d0; 
595  if (d == 0)
596  {
597    if (substitute > 1)
598      return N (reverseSubst (gcdcAcB, substitute));
599    else 
600      return N(gcdcAcB);
601  }
602
603  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
604  CanonicalForm newtonPoly;
605
606  newtonPoly= 1;
607  m= gcdlcAlcB;
608  G_m= 0;
609  H= 0;
610  bool fail= false;
611  topLevel= false;
612  bool inextension= false;
613  Variable V_buf= alpha;
614  CanonicalForm prim_elem, im_prim_elem;
615  CFList source, dest;
616  do 
617  {
618    random_element= randomElement (m, V_buf, l, fail);
619    if (fail) 
620    {
621      source= CFList();
622      dest= CFList();
623      int num_vars= tmin (getNumVars(A), getNumVars(B));
624      num_vars--;
625
626      choose_extension (d, num_vars, V_buf);
627      bool prim_fail= false;
628      Variable V_buf2;
629      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
630
631      ASSERT (!prim_fail, "failure in integer factorizer");
632      if (prim_fail)
633        ; //ERROR
634      else
635        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
636
637      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
638      DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
639      inextension= true;
640      for (CFListIterator i= l; i.hasItem(); i++) 
641        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem, 
642                             im_prim_elem, source, dest);
643      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
644      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
645      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, 
646                          source, dest);
647      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
648      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
649      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, 
650                         source, dest);
651
652      fail= false;
653      random_element= randomElement (m, V_buf, l, fail );
654      DEBOUTLN (cerr, "fail= " << fail);
655      CFList list;
656      TIMING_START (gcd_recursion);
657      G_random_element= 
658      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
659                        list, topLevel);
660      TIMING_END_AND_PRINT (gcd_recursion, 
661                            "time for recursive call: ");
662      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
663    }
664    else 
665    {
666      CFList list;
667      TIMING_START (gcd_recursion);
668      G_random_element= 
669      GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf,
670                        list, topLevel);
671      TIMING_END_AND_PRINT (gcd_recursion, 
672                            "time for recursive call: ");
673      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
674    }
675
676    d0= totaldegree (G_random_element, Variable(2), 
677                     Variable (G_random_element.level()));
678    if (d0 == 0)
679    {
680      if (substitute > 1)
681        return N (reverseSubst (gcdcAcB, substitute));
682      else 
683        return N(gcdcAcB);
684    }
685    if (d0 >  d) 
686    {
687      if (!find (l, random_element))
688        l.append (random_element);
689      continue;
690    }
691
692    G_random_element= 
693    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
694    * G_random_element;
695
696    d0= totaldegree (G_random_element, Variable(2), 
697                     Variable(G_random_element.level()));
698    if (d0 <  d) 
699    {
700      m= gcdlcAlcB;
701      newtonPoly= 1;
702      G_m= 0;
703      d= d0;
704    }
705
706    TIMING_START (newton_interpolation);
707    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
708    TIMING_END_AND_PRINT (newton_interpolation, 
709                          "time for newton interpolation: ");
710
711    //termination test
712    if (uni_lcoeff (H) == gcdlcAlcB) 
713    {
714      cH= uni_content (H);
715      ppH= H/cH;
716      if (inextension) 
717      {
718        CFList u, v; 
719        //maybe it's better to test if ppH is an element of F(\alpha) before
720        //mapping down
721        DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
722        ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v); 
723        ppH /= Lc(ppH);
724        DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
725        if (fdivides (ppH, A) && fdivides (ppH, B))         
726        {
727          if (substitute > 1)
728          {
729            ppH= reverseSubst (ppH, substitute);
730            gcdcAcB= reverseSubst (gcdcAcB, substitute);
731          }
732          return N(gcdcAcB*ppH);
733        }
734      }
735      else if (fdivides (ppH, A) && fdivides (ppH, B)) 
736      {
737        if (substitute > 1)
738        {
739          ppH= reverseSubst (ppH, substitute);
740          gcdcAcB= reverseSubst (gcdcAcB, substitute);
741        }
742        return N(gcdcAcB*ppH);
743      }
744    }
745
746    G_m= H;
747    newtonPoly= newtonPoly*(x - random_element);
748    m= m*(x - random_element);
749    if (!find (l, random_element))
750      l.append (random_element);
751  } while (1);
752}
753
754/// compute a random element a of GF, s.t. F(a) \f$ \neq 0 \f$ , F is a
755/// univariate polynomial, returns fail if there are no field elements left
756/// which have not been used before
757static inline
758CanonicalForm
759GFRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
760{
761  fail= false;
762  Variable x= F.mvar();
763  GFRandom genGF;
764  CanonicalForm random;
765  int p= getCharacteristic();
766  int d= getGFDegree();
767  int bound= ipower (p, d);
768  do 
769  {
770    if (list.length() == bound)
771    {
772      fail= true;
773      break;
774    }
775    if (list.length() < 1)
776      random= 0;
777    else 
778    {
779      random= genGF.generate();
780      while (find (list, random))
781        random= genGF.generate();
782    }
783    if (F (random, x) == 0) 
784    {
785      list.append (random);
786      continue;
787    }
788  } while (find (list, random));
789  return random;
790}
791
792/// GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for
793/// Computer Algebra" by Geddes, Czapor, Labahn
794/// Usually this algorithm will be faster than GCD_Fp_extension since GF has
795/// faster field arithmetics, however it might fail if the input is large since
796/// the size of the base field is bounded by 2^16, output is monic
797CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
798        CFList& l, bool& topLevel) 
799{ 
800  CanonicalForm A= F;
801  CanonicalForm B= G;
802  if (F.isZero() && degree(G) > 0) return G/Lc(G);
803  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
804  else if (F.isZero() && G.isZero()) return F.genOne();
805  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
806  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
807  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
808  if (F == G) return F/Lc(F);
809 
810  CFMap M,N;
811  int best_level= myCompress (A, B, M, N, topLevel);
812
813  if (best_level == 0) return B.genOne();
814
815  A= M(A);
816  B= M(B);
817
818  Variable x= Variable(1);
819
820  //univariate case
821  if (A.isUnivariate() && B.isUnivariate()) 
822    return N (gcd (A, B)); 
823
824  int substitute= substituteCheck (A, B);
825 
826  if (substitute > 1)
827    subst (A, B, A, B, substitute);
828
829  CanonicalForm cA, cB;    // content of A and B
830  CanonicalForm ppA, ppB;    // primitive part of A and B
831  CanonicalForm gcdcAcB;
832
833  if (topLevel)
834  {
835    if (best_level <= 2)
836      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
837    else 
838      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 
839  }
840  else
841  {
842    cA = uni_content (A);
843    cB = uni_content (B); 
844    gcdcAcB= gcd (cA, cB);
845    ppA= A/cA;
846    ppB= B/cB;
847  }
848
849  CanonicalForm lcA, lcB;  // leading coefficients of A and B
850  CanonicalForm gcdlcAlcB; 
851
852  lcA= uni_lcoeff (ppA);
853  lcB= uni_lcoeff (ppB);
854
855  if (fdivides (lcA, lcB)) 
856  { 
857    if (fdivides (A, B))
858      return F; 
859  }   
860  if (fdivides (lcB, lcA)) 
861  { 
862    if (fdivides (B, A)) 
863      return G;
864  }
865
866  gcdlcAlcB= gcd (lcA, lcB);
867
868  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
869  if (d == 0)
870  {
871    if (substitute > 1)
872      return N (reverseSubst (gcdcAcB, substitute));
873    else 
874      return N(gcdcAcB);
875  }
876  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
877  if (d0 < d) 
878    d= d0; 
879  if (d == 0)
880  {
881    if (substitute > 1)
882      return N (reverseSubst (gcdcAcB, substitute));
883    else 
884      return N(gcdcAcB);
885  }
886
887  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
888  CanonicalForm newtonPoly;
889
890  newtonPoly= 1;
891  m= gcdlcAlcB;
892  G_m= 0;
893  H= 0;
894  bool fail= false;
895  topLevel= false;
896  bool inextension= false;
897  int p;
898  int k= getGFDegree();
899  int kk;
900  int expon;
901  char gf_name_buf= gf_name;
902  do 
903  {
904    random_element= GFRandomElement (m, l, fail);
905    if (fail) 
906    { 
907      int num_vars= tmin (getNumVars(A), getNumVars(B));
908      num_vars--;
909      p= getCharacteristic();
910      expon= (int) floor ((log ((double) ipower (d + 1, num_vars))/log ((double) p)));
911      if (expon < 2)
912        expon= 2;
913      kk= getGFDegree(); 
914      if (ipower (p, kk*expon) < (1 << 16)) 
915        setCharacteristic (p, kk*(int)expon, 'b');
916      else 
917      {
918        expon= (int) floor((log ((double)((1<<16) - 1)))/(log((double)p)*kk));
919        ASSERT (expon >= 2, "not enough points in GCD_GF");
920        setCharacteristic (p, (int)(kk*expon), 'b');
921      }
922      inextension= true;
923      fail= false;
924      for (CFListIterator i= l; i.hasItem(); i++) 
925        i.getItem()= GFMapUp (i.getItem(), kk);
926      m= GFMapUp (m, kk);
927      G_m= GFMapUp (G_m, kk);
928      newtonPoly= GFMapUp (newtonPoly, kk);
929      ppA= GFMapUp (ppA, kk);
930      ppB= GFMapUp (ppB, kk);
931      gcdlcAlcB= GFMapUp (gcdlcAlcB, kk);
932      random_element= GFRandomElement (m, l, fail);
933      DEBOUTLN (cerr, "fail= " << fail);
934      CFList list;
935      TIMING_START (gcd_recursion);
936      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 
937                                list, topLevel);
938      TIMING_END_AND_PRINT (gcd_recursion, 
939                            "time for recursive call: ");
940      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
941    }
942    else 
943    {
944      CFList list;
945      TIMING_START (gcd_recursion);
946      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x), 
947                                list, topLevel);
948      TIMING_END_AND_PRINT (gcd_recursion, 
949                            "time for recursive call: ");
950      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
951    }
952
953    d0= totaldegree (G_random_element, Variable(2), 
954                     Variable (G_random_element.level()));
955    if (d0 == 0) 
956    {
957      if (inextension) 
958      {
959        setCharacteristic (p, k, gf_name_buf);
960        {
961          if (substitute > 1)
962            return N (reverseSubst (gcdcAcB, substitute));
963          else
964            return N(gcdcAcB); 
965        } 
966      }
967      else
968      {
969        if (substitute > 1)
970          return N (reverseSubst (gcdcAcB, substitute));
971        else
972          return N(gcdcAcB); 
973      }
974    } 
975    if (d0 >  d) 
976    {
977      if (!find (l, random_element))
978        l.append (random_element);
979      continue;
980    }
981
982    G_random_element= 
983    (gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) *
984     G_random_element;
985    d0= totaldegree (G_random_element, Variable(2), 
986                     Variable (G_random_element.level()));
987
988    if (d0 < d) 
989    {
990      m= gcdlcAlcB;
991      newtonPoly= 1;
992      G_m= 0;
993      d= d0;
994    }
995
996    TIMING_START (newton_interpolation);
997    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
998    TIMING_END_AND_PRINT (newton_interpolation, "time for newton interpolation: ");
999
1000    //termination test
1001    if (uni_lcoeff (H) == gcdlcAlcB) 
1002    {
1003      cH= uni_content (H);
1004      ppH= H/cH;
1005      if (inextension) 
1006      {
1007        if (fdivides(ppH, GFMapUp(A, k)) && fdivides(ppH, GFMapUp(B,k))) 
1008        {
1009          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
1010          ppH= GFMapDown (ppH, k);
1011          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
1012          if (substitute > 1)
1013          {
1014            ppH= reverseSubst (ppH, substitute);
1015            gcdcAcB= reverseSubst (gcdcAcB, substitute);
1016          }
1017          setCharacteristic (p, k, gf_name_buf);
1018          return N(gcdcAcB*ppH);
1019        }
1020      }
1021      else 
1022      {
1023        if (fdivides (ppH, A) && fdivides (ppH, B)) 
1024        {
1025          if (substitute > 1)
1026          {
1027            ppH= reverseSubst (ppH, substitute);
1028            gcdcAcB= reverseSubst (gcdcAcB, substitute);
1029          }
1030          return N(gcdcAcB*ppH);
1031        }
1032      }
1033    }
1034
1035    G_m= H;
1036    newtonPoly= newtonPoly*(x - random_element);
1037    m= m*(x - random_element);
1038    if (!find (l, random_element))
1039      l.append (random_element);
1040  } while (1);
1041}
1042
1043/// F is assumed to be an univariate polynomial in x,
1044/// computes a random monic irreducible univariate polynomial of random
1045/// degree < i in x which does not divide F
1046CanonicalForm
1047randomIrredpoly (int i, const Variable & x) 
1048{
1049  int p= getCharacteristic();
1050  ZZ NTLp= to_ZZ (p);
1051  ZZ_p::init (NTLp);
1052  ZZ_pX NTLirredpoly; 
1053  CanonicalForm CFirredpoly;
1054  BuildIrred (NTLirredpoly, i + 1);
1055  CFirredpoly= convertNTLZZpX2CF (NTLirredpoly, x);
1056  return CFirredpoly;
1057} 
1058
1059static inline
1060CanonicalForm
1061FpRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
1062{
1063  fail= false;
1064  Variable x= F.mvar();
1065  FFRandom genFF;
1066  CanonicalForm random;
1067  int p= getCharacteristic();
1068  int bound= p;
1069  do 
1070  {
1071    if (list.length() == bound)
1072    {
1073      fail= true;
1074      break;
1075    }
1076    if (list.length() < 1)
1077      random= 0;
1078    else 
1079    {
1080      random= genFF.generate();
1081      while (find (list, random))
1082        random= genFF.generate();
1083    }
1084    if (F (random, x) == 0) 
1085    {
1086      list.append (random);
1087      continue;
1088    }
1089  } while (find (list, random));
1090  return random;
1091}
1092
1093CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G, 
1094                           bool& topLevel, CFList& l) 
1095{
1096  CanonicalForm A= F;
1097  CanonicalForm B= G;
1098  if (F.isZero() && degree(G) > 0) return G/Lc(G);
1099  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
1100  else if (F.isZero() && G.isZero()) return F.genOne();
1101  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
1102  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
1103  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
1104  if (F == G) return F/Lc(F);
1105
1106  CFMap M,N;
1107  int best_level= myCompress (A, B, M, N, topLevel);
1108
1109  if (best_level == 0) return B.genOne();
1110
1111  A= M(A);
1112  B= M(B);
1113
1114  Variable x= Variable (1);
1115
1116  //univariate case
1117  if (A.isUnivariate() && B.isUnivariate()) 
1118    return N (gcd (A, B));
1119
1120  int substitute= substituteCheck (A, B);
1121 
1122  if (substitute > 1)
1123    subst (A, B, A, B, substitute);
1124
1125  CanonicalForm cA, cB;    // content of A and B
1126  CanonicalForm ppA, ppB;    // primitive part of A and B
1127  CanonicalForm gcdcAcB;
1128
1129  if (topLevel)
1130  {
1131    if (best_level <= 2)
1132      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
1133    else 
1134      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); 
1135  }
1136  else
1137  {
1138    cA = uni_content (A);
1139    cB = uni_content (B); 
1140    gcdcAcB= gcd (cA, cB);
1141    ppA= A/cA;
1142    ppB= B/cB;
1143  }
1144
1145  CanonicalForm lcA, lcB;  // leading coefficients of A and B
1146  CanonicalForm gcdlcAlcB; 
1147  lcA= uni_lcoeff (ppA);
1148  lcB= uni_lcoeff (ppB);
1149
1150  if (fdivides (lcA, lcB)) 
1151  { 
1152    if (fdivides (A, B))
1153      return F/Lc(F);
1154  }   
1155  if (fdivides (lcB, lcA))
1156  { 
1157    if (fdivides (B, A)) 
1158      return G/Lc(G);
1159  }
1160
1161  gcdlcAlcB= gcd (lcA, lcB); 
1162 
1163  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
1164  int d0;
1165
1166  if (d == 0)
1167  {
1168    if (substitute > 1)
1169      return N (reverseSubst (gcdcAcB, substitute));
1170    else 
1171      return N(gcdcAcB);
1172  }
1173  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
1174
1175  if (d0 < d) 
1176    d= d0;
1177
1178  if (d == 0) 
1179  {
1180    if (substitute > 1)
1181      return N (reverseSubst (gcdcAcB, substitute));
1182    else 
1183      return N(gcdcAcB);
1184  }
1185
1186  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
1187  CanonicalForm newtonPoly= 1;
1188  m= gcdlcAlcB;
1189  H= 0;
1190  G_m= 0;
1191  Variable alpha, V_buf;
1192  bool fail= false;
1193  bool inextension= false;
1194  bool inextensionextension= false;
1195  topLevel= false;
1196  CFList source, dest;
1197  do 
1198  {
1199    if (inextension)
1200      random_element= randomElement (m, alpha, l, fail);
1201    else
1202      random_element= FpRandomElement (m, l, fail);
1203
1204    if (!fail && !inextension)
1205    {
1206      CFList list;
1207      TIMING_START (gcd_recursion);
1208      G_random_element= 
1209      GCD_small_p (ppA (random_element,x), ppB (random_element,x), topLevel, 
1210      list);
1211      TIMING_END_AND_PRINT (gcd_recursion, 
1212                            "time for recursive call: ");
1213      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1214    }
1215    else if (!fail && inextension)
1216    {
1217      CFList list; 
1218      TIMING_START (gcd_recursion);
1219      G_random_element= 
1220      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 
1221                        list, topLevel);
1222      TIMING_END_AND_PRINT (gcd_recursion, 
1223                            "time for recursive call: ");
1224      DEBOUTLN (cerr, "G_random_element= " << G_random_element);   
1225    }
1226    else if (fail && !inextension)
1227    {
1228      source= CFList();
1229      dest= CFList();
1230      CFList list;
1231      CanonicalForm mipo;
1232      int deg= 2;
1233      do {
1234        mipo= randomIrredpoly (deg, x); 
1235        alpha= rootOf (mipo);
1236        inextension= true;
1237        fail= false;
1238        random_element= randomElement (m, alpha, l, fail); 
1239        deg++;
1240      } while (fail); 
1241      list= CFList();
1242      TIMING_START (gcd_recursion);
1243      G_random_element= 
1244      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha, 
1245                        list, topLevel);
1246      TIMING_END_AND_PRINT (gcd_recursion, 
1247                            "time for recursive call: ");
1248      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1249    }
1250    else if (fail && inextension)
1251    {
1252      source= CFList();
1253      dest= CFList();
1254      int num_vars= tmin (getNumVars(A), getNumVars(B));
1255      num_vars--;
1256      V_buf= alpha;
1257      choose_extension (d, num_vars, V_buf);
1258      bool prim_fail= false;
1259      Variable V_buf2;
1260      CanonicalForm prim_elem, im_prim_elem;
1261      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
1262       
1263      ASSERT (!prim_fail, "failure in integer factorizer");
1264      if (prim_fail)
1265        ; //ERROR
1266      else
1267        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
1268
1269      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
1270      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
1271
1272      inextensionextension= true;
1273      for (CFListIterator i= l; i.hasItem(); i++) 
1274        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem, 
1275                             im_prim_elem, source, dest);
1276      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1277      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1278      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem, 
1279                          source, dest);
1280      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1281      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1282      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem, 
1283                         source, dest);
1284      fail= false;
1285      random_element= randomElement (m, V_buf, l, fail );
1286      DEBOUTLN (cerr, "fail= " << fail);
1287      CFList list;
1288      TIMING_START (gcd_recursion);
1289      G_random_element= 
1290      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
1291                        list, topLevel);
1292      TIMING_END_AND_PRINT (gcd_recursion, 
1293                            "time for recursive call: ");
1294      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1295      alpha= V_buf;
1296    } 
1297
1298    d0= totaldegree (G_random_element, Variable(2), 
1299                     Variable (G_random_element.level()));
1300
1301    if (d0 == 0)
1302    {
1303      if (substitute > 1)
1304        return N (reverseSubst (gcdcAcB, substitute));
1305      else
1306        return N(gcdcAcB); 
1307    }
1308    if (d0 >  d) 
1309    { 
1310      if (!find (l, random_element))
1311        l.append (random_element);
1312      continue;
1313    }
1314
1315    G_random_element= (gcdlcAlcB(random_element,x)/uni_lcoeff(G_random_element))
1316                       *G_random_element; 
1317
1318   
1319    d0= totaldegree (G_random_element, Variable(2), 
1320                     Variable(G_random_element.level()));
1321
1322    if (d0 <  d) 
1323    {
1324      m= gcdlcAlcB;
1325      newtonPoly= 1;
1326      G_m= 0;
1327      d= d0;
1328    }
1329     
1330    TIMING_START (newton_interpolation);
1331    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1332    TIMING_END_AND_PRINT (newton_interpolation, 
1333                          "time for newton_interpolation: ");
1334
1335    //termination test
1336    if (uni_lcoeff (H) == gcdlcAlcB) 
1337    {
1338      cH= uni_content (H);
1339      ppH= H/cH;
1340      ppH /= Lc (ppH);
1341      DEBOUTLN (cerr, "ppH= " << ppH);
1342      if (fdivides (ppH, A) && fdivides (ppH, B)) 
1343      {
1344        if (substitute > 1)
1345        {
1346          ppH= reverseSubst (ppH, substitute);
1347          gcdcAcB= reverseSubst (gcdcAcB, substitute);
1348        }
1349        return N(gcdcAcB*ppH);
1350      }
1351    }
1352
1353    G_m= H;
1354    newtonPoly= newtonPoly*(x - random_element);
1355    m= m*(x - random_element);
1356    if (!find (l, random_element))
1357      l.append (random_element);
1358  } while (1);
1359}
1360
1361#endif
Note: See TracBrowser for help on using the repository browser.