source: git/factory/facFqFactorize.cc @ 883ea0b

spielwiese
Last change on this file since 883ea0b was 883ea0b, checked in by Hans Schoenemann <hannes@…>, 14 years ago
ceil log -> ilog2 git-svn-id: file:///usr/local/Singular/svn/trunk@12926 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 53.8 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facFqFactorize.cc
5 *
6 * This file implements functions for factoring a multivariate polynomial over
7 * a finite field.
8 *
9 * ABSTRACT: "Efficient Multivariate Factorization over Finite Fields" by
10 * L. Bernardin & M. Monagon.
11 *
12 * @author Martin Lee
13 *
14 * @internal @version \$Id$
15 *
16 **/
17/*****************************************************************************/
18
19#include <config.h>
20
21#include "assert.h"
22#include "debug.h"
23#include "timing.h"
24
25#include "facFqFactorizeUtil.h"
26#include "facFqFactorize.h"
27#include "cf_random.h"
28#include "facHensel.h"
29#include "cf_gcd_smallp.h"
30#include "cf_map_ext.h"
31#include "algext.h"
32
33#ifdef HAVE_NTL
34#include <NTL/ZZ_pEX.h>
35#include "NTLconvert.h"
36
37TIMING_DEFINE_PRINT(fac_bi_factorizer);
38TIMING_DEFINE_PRINT(fac_hensel_lift);
39TIMING_DEFINE_PRINT(fac_factor_recombination);
40
41static inline
42CanonicalForm
43listGCD (const CFList& L);
44
45static inline
46CanonicalForm
47myContent (const CanonicalForm& F)
48{
49  CanonicalForm G= swapvar (F, F.mvar(), Variable (1));
50  CFList L;
51  Variable alpha;
52  bool algExt= hasFirstAlgVar (G, alpha);
53  for (CFIterator i= G; i.hasTerms(); i++)
54    L.append (i.coeff());
55  if (L.length() == 2)
56  {
57    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
58    if (GF)
59    {
60      return swapvar (GCD_GF (L.getFirst(), L.getLast()), F.mvar(),
61                      Variable (1));
62    }
63    else if (!GF && algExt)
64    {
65      return swapvar (GCD_Fp_extension (L.getFirst(), L.getLast(), alpha),
66                                       F.mvar(), Variable (1));
67    }
68    else
69      return swapvar (GCD_small_p (L.getFirst(), L.getLast()), F.mvar(),
70                     Variable (1));
71  }
72  if (L.length() == 1)
73  {
74    return LC (F, Variable (1));
75  }
76  return swapvar (listGCD (L), F.mvar(), Variable (1));
77}
78
79static inline
80CanonicalForm
81listGCD (const CFList& L)
82{
83  if (L.length() == 1)
84    return L.getFirst();
85  Variable alpha;
86  if (L.length() == 2)
87  {
88    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
89    if (GF)
90    {
91      return GCD_GF (L.getFirst(), L.getLast());
92    }
93    else if (!GF && (hasFirstAlgVar (L.getFirst(), alpha) ||
94                     hasFirstAlgVar (L.getLast(), alpha)))
95    {
96      return GCD_Fp_extension (L.getFirst(), L.getLast(), alpha);
97    }
98    else
99      return GCD_small_p (L.getFirst(), L.getLast());
100  }
101  else
102  {
103    CFList lHi, lLo;
104    CanonicalForm resultHi, resultLo;
105    int length= L.length()/2;
106    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
107    int j= 0;
108    for (CFListIterator i= L; j < length; i++, j++)
109      lHi.append (i.getItem());
110    lLo= Difference (L, lHi);
111    resultHi= listGCD (lHi);
112    resultLo= listGCD (lLo);
113    if (resultHi.isOne() || resultLo.isOne())
114      return 1;
115    if (GF)
116    {
117      return GCD_GF (resultHi, resultLo);
118    }
119    else if (!GF && (hasFirstAlgVar (resultHi, alpha) ||
120                     hasFirstAlgVar (resultLo, alpha)))
121    {
122      return GCD_Fp_extension (resultHi, resultLo, alpha);
123    }
124    else
125      return GCD_small_p (resultHi, resultLo);
126  }
127}
128
129static inline
130CanonicalForm
131myContent (const CanonicalForm& F, const Variable& x)
132{
133  if (degree (F, x) <= 0)
134    return 1;
135  CanonicalForm G= F;
136  bool swap= false;
137  if (x != F.mvar())
138  {
139    swap= true;
140    G= swapvar (F, x, F.mvar());
141  }
142  CFList L;
143  Variable alpha;
144  bool algExt= hasFirstAlgVar (G, alpha);
145  for (CFIterator i= G; i.hasTerms(); i++)
146    L.append (i.coeff());
147  if (L.length() == 2)
148  {
149    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
150    if (GF)
151    {
152      if (swap)
153        return swapvar(GCD_GF (L.getFirst(), L.getLast()), F.mvar(), x);
154      else
155        return GCD_GF (L.getFirst(), L.getLast());
156    }
157    else if (!GF && algExt)
158    {
159      if (swap)
160        return swapvar(GCD_Fp_extension (L.getFirst(), L.getLast(), alpha),
161                                         F.mvar(), x);
162      else
163        return GCD_Fp_extension (L.getFirst(), L.getLast(), alpha);
164    }
165    else
166    {
167      if (swap)
168        return swapvar (GCD_small_p (L.getFirst(), L.getLast()), F.mvar(),
169                        x);
170      else
171        return GCD_small_p (L.getFirst(), L.getLast());
172    }
173  }
174  if (L.length() == 1)
175  {
176    return LC (F, x);
177  }
178  if (swap)
179    return swapvar (listGCD (L), F.mvar(), x);
180  else
181    return listGCD (L);
182}
183
184static inline
185CanonicalForm
186myLcm (const CanonicalForm& F, const CanonicalForm& G)
187{
188  if ( F.isZero() || G.isZero() )
189    return 0;
190  else
191  {
192    Variable alpha;
193    bool GF= (CFFactory::gettype() == GaloisFieldDomain);
194    if (GF)
195      return (F/GCD_GF (F, G))*G;
196    else if (!GF && (hasFirstAlgVar (F, alpha) ||
197                     hasFirstAlgVar (G, alpha)))
198      return (F/GCD_Fp_extension (F, G, alpha))*G;
199    else
200      return (F/GCD_small_p (F, G))*G;
201  }
202}
203
204
205static inline
206CanonicalForm myCompress (const CanonicalForm& F, CFMap& N)
207{
208  int n= F.level();
209  int * degsf= new int [n + 1];
210  int ** swap;
211  swap= new int* [n + 1];
212  for (int i= 0; i <= n; i++)
213  {
214    degsf[i]= 0;
215    swap [i]= new int [2];
216    swap [i] [0]= 0;
217    swap [i] [1]= 0;
218  }
219  int i= 1;
220  n= 1;
221  degsf= degrees (F, degsf);
222
223  CanonicalForm result= F;
224  while ( i <= F.level() )
225  {
226    while( degsf[i] == 0 ) i++;
227    swap[n][0]= i;
228    swap[n][1]= degsf[i];
229    if (i != n)
230      result= swapvar (result, Variable (n), Variable(i));
231    n++; i++;
232  }
233
234  int buf1, buf2;
235  n--;
236
237  for (i= 1; i < n; i++)
238  {
239    for (int j= 1; j < n - i + 1; j++)
240    {
241      if (swap[j][1] < swap[j + 1][1])
242      {
243        buf1= swap [j + 1] [0];
244        buf2= swap [j + 1] [1];
245        swap[j + 1] [0]= swap[j] [0];
246        swap[j + 1] [1]= swap[j] [1];
247        swap[j][0]= buf1;
248        swap[j][1]= buf2;
249        result= swapvar (result, Variable (j + 1), Variable (j));
250      }
251    }
252  }
253
254  for (i= n; i > 0; i--)
255  {
256    if (i != swap[i] [0])
257      N.newpair (Variable (i), Variable (swap[i] [0]));
258  }
259
260  for (i= 0; i <= n; i++)
261    delete [] swap[i];
262  delete [] swap;
263
264  delete [] degsf;
265
266  return result;
267}
268
269CFList
270monicFactorRecombi (const CFList& factors,const CanonicalForm& F, const
271                    CanonicalForm& M, const DegreePattern& degs)
272{
273  if (degs.getLength() == 1)
274    return CFList(F);
275
276  CFList T, S;
277
278  T= factors;
279
280  int s= 1;
281  CFList result;
282  CanonicalForm buf= F;
283  CanonicalForm LCBuf= LC (buf, Variable (1));
284  CanonicalForm g, gg;
285  int v [T.length()];
286  for (int i= 0; i < T.length(); i++)
287    v[i]= 0;
288  bool noSubset= false;
289  CFArray TT;
290  int subsetDeg;
291  DegreePattern bufDegs1= degs, bufDegs2;
292  TT= copy (factors);
293  while (T.length() >= 2*s)
294  {
295    while (noSubset == false)
296    {
297      if (T.length() == s)
298      {
299        result.append (prodMod (T, M));
300        return result;
301      }
302      S= subset (v, s, TT, noSubset);
303      if (noSubset) break;
304      subsetDeg= subsetDegree (S);
305      if (!degs.find (subsetDeg))
306        continue;
307      else
308      {
309        g= prodMod0 (S, M);
310        gg= mod (g*LCBuf, M); //univariate polys
311        gg /= content (gg);
312        if (fdivides (LC (gg), LCBuf))
313        {
314          g= prodMod (S, M);
315          gg= mulMod2 (LCBuf, g, M);
316          gg /= content (gg, Variable (1));
317          if (fdivides (gg, buf))
318          {
319            result.append (g);
320            buf /= gg;
321            LCBuf= LC (buf, Variable(1));
322            T= Difference (T, S);
323            // compute new possible degree pattern
324            bufDegs2= DegreePattern (T);
325            bufDegs1.intersect (bufDegs2);
326            bufDegs1.refine ();
327            if (T.length() < 2*s || T.length() == s ||
328                bufDegs1.getLength() == 1)
329            {
330              result.append (prodMod (T, M));
331              return result;
332            }
333            TT= copy (T);
334            indexUpdate (v, s, T.length(), noSubset);
335            if (noSubset) break;
336          }
337        }
338      }
339    }
340    s++;
341    if (T.length() < 2*s || T.length() == s)
342    {
343      result.append (prodMod (T, M));
344      return result;
345    }
346    int v [T.length()];
347    for (int i= 0; i < T.length(); i++)
348      v[i]= 0;
349    noSubset= false;
350  }
351  if (T.length() < 2*s)
352    result.append (prodMod (T, M));
353
354  return result;
355}
356
357CFList
358earlyMonicFactorDetect (CanonicalForm& F, CFList& factors,
359                        int& adaptedLiftBound, DegreePattern& degs,
360                        bool& success, int deg, const int bound)
361{
362  DegreePattern bufDegs1= degs;
363  DegreePattern bufDegs2;
364  CFList result;
365  CFList T= factors;
366  CanonicalForm buf= F;
367  CanonicalForm LCBuf= LC (buf, Variable (1));
368  CanonicalForm g, gg;
369  CanonicalForm M= power (F.mvar(), deg);
370  int d= bound;
371  int e= 0;
372  adaptedLiftBound= 0;
373  for (CFListIterator i= factors; i.hasItem(); i++)
374  {
375    if (!bufDegs1.find (degree (i.getItem(), 1)))
376      continue;
377    else
378    {
379      gg= i.getItem() (0, 1);
380      gg *= LCBuf;
381      gg= mod (gg, M);
382      if (fdivides (LC (gg), LCBuf))
383      {
384        g= i.getItem();
385        gg= mulMod2 (g, LCBuf, M);
386        gg /= content (gg, Variable (1));
387        if (fdivides (gg, buf))
388        {
389          result.append (g);
390          CanonicalForm bufbuf= buf;
391          buf /= gg;
392          d -= degree (gg) + degree (LC (gg, 1));
393          e= tmax (e, degree (gg) + degree (LC (gg, 1)));
394          LCBuf= LC (buf, Variable (1));
395          T= Difference (T, CFList (i.getItem()));
396
397          // compute new possible degree pattern
398          bufDegs2= DegreePattern (T);
399          bufDegs1.intersect (bufDegs2);
400          bufDegs1.refine ();
401          if (bufDegs1.getLength() <= 1)
402          {
403            result.append (prodMod (T, M));
404            break;
405          }
406        }
407      }
408    }
409  }
410  adaptedLiftBound= d;
411
412  if (adaptedLiftBound < deg)
413  {
414    factors= T;
415    degs= bufDegs1;
416    if (adaptedLiftBound < degree (F) + 1)
417    {
418      if (d == 1)
419      {
420        if (e + 1 > deg)
421        {
422          success= false;
423          adaptedLiftBound= deg;
424        }
425        else
426        {
427          F= buf;
428          success= true;
429          if (e + 1 < degree (F) + 1)
430            adaptedLiftBound= deg;
431          else
432            adaptedLiftBound= e + 1;
433        }
434      }
435      else
436      {
437        success= true;
438        F= buf;
439        adaptedLiftBound= deg;
440      }
441    }
442    else
443    {
444      F= buf;
445      success= true;
446    }
447  }
448  if (bufDegs1.getLength() <= 1)
449  {
450    degs= bufDegs1;
451    if (!success)
452      adaptedLiftBound= deg;
453  }
454
455  return result;
456}
457
458CFList biFactorizer (const CanonicalForm& F, const Variable& alpha,
459                     CanonicalForm& bivarEval, int& liftBound)
460{
461  CanonicalForm A= F;
462  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
463
464  if (A.isUnivariate())
465    return uniFactorizer (F, alpha, GF);
466
467
468  CFMap N;
469  A= compress (A, N);
470  Variable y= A.mvar();
471  A /= Lc(A);
472
473  if (y.level() > 2) return CFList (F);
474  Variable x= Variable (1);
475
476  //remove content and factorize content
477  CanonicalForm contentAy= content (A, y);
478  CanonicalForm contentAx= content (A, x);
479  A= A/(contentAy*contentAx);
480  CFList contentAyFactors= uniFactorizer (contentAy, alpha, GF);
481
482  //trivial case
483  CFList factors;
484  if (A.inCoeffDomain())
485  {
486    append (factors, contentAyFactors);
487    decompress (factors, N);
488    bivarEval= N (y - bivarEval);
489    return factors;
490  }
491  else if (A.isUnivariate())
492  {
493    if (A.mvar() == x)
494      factors= uniFactorizer (A, alpha, GF);
495    append (factors, contentAyFactors);
496    decompress (factors, N);
497    bivarEval= N (y - bivarEval);
498    return factors;
499  }
500  bool fail;
501  CanonicalForm Aeval, evaluation, bufAeval, bufEvaluation, buf;
502  CFList uniFactors, list, bufUniFactors;
503  DegreePattern degs;
504  DegreePattern bufDegs;
505
506  int factorNums= 5;
507  if (factorNums < (int) ilog2 (totaldegree (A)))
508    factorNums= (int) ilog2 (totaldegree (A));
509  for (int i= 0; i < factorNums; i++)
510  {
511    if (i == 0)
512      Aeval= A (bivarEval, y);
513    else if (i > 0 && contentAx.inCoeffDomain())
514    {
515      Aeval= A;
516      evaluation= evalPoint (A, Aeval, alpha, list, GF, fail);
517    }
518
519    if (fail && (i != 0))
520      break;
521
522    // univariate factorization
523    uniFactors= uniFactorizer (Aeval, alpha, GF);
524
525    if (uniFactors.length() == 1)
526    {
527      append (factors, contentAyFactors);
528      if (contentAyFactors.isEmpty())
529        factors.append (F/lc(F));
530      else
531      {
532        A= A (y - bivarEval, y);
533        A= A/lc (A);
534        if (!LC(A, 1).isOne())
535        {
536          CanonicalForm s,t;
537          (void) extgcd (LC (A, 1), power (y, liftBound), s, t);
538          A *= s;
539          A= mod (A, power (y, liftBound));
540        }
541        factors.append (A);
542      }
543      decompress (factors, N);
544      bivarEval= N (bivarEval);
545      return factors;
546    }
547
548    // degree analysis
549    degs= DegreePattern (uniFactors);
550
551    if (i == 0)
552    {
553      bufAeval= Aeval;
554      bufEvaluation= bivarEval;
555      bufUniFactors= uniFactors;
556      bufDegs= degs;
557      if (!contentAx.inCoeffDomain())
558        break;
559    }
560    else
561    {
562      bufDegs.intersect (degs);
563      if (uniFactors.length() < bufUniFactors.length())
564      {
565        bufUniFactors= uniFactors;
566        bufAeval= Aeval;
567        bufEvaluation= evaluation;
568      }
569    }
570
571    if (bufDegs.getLength() == 1)
572    {
573      append (factors, contentAyFactors);
574      if (contentAyFactors.isEmpty())
575        factors.append (F/lc(F));
576      else
577      {
578        A= A (y - bivarEval, y);
579        A= A/lc (A);
580        if (!LC(A, 1).isOne())
581        {
582          CanonicalForm s,t;
583          (void) extgcd (LC (A, 1), power (y, liftBound), s, t);
584          A *= s;
585          A= mod (A, power (y, liftBound));
586        }
587        factors.append (A);
588      }
589      decompress (factors, N);
590      bivarEval= N (bivarEval);
591      return factors;
592    }
593    list.append (evaluation);
594  }
595
596  bivarEval= y - bufEvaluation;
597  A= A (y + bufEvaluation, y);
598  bufUniFactors.insert (LC (A, x));
599
600  // Hensel lifting
601  CFList diophant;
602  CFMatrix M= CFMatrix (liftBound, bufUniFactors.length() - 1);
603  CFArray Pi;
604  bool earlySuccess= false;
605  int newLiftBound= 0;
606  CFList earlyFactors;
607  int smallFactorDeg= 11; //tunable parameter
608  if (smallFactorDeg >= liftBound)
609    henselLift12 (A, bufUniFactors, liftBound, Pi, diophant, M);
610  else if (smallFactorDeg >= degree (A, y) + 1)
611  {
612    henselLift12 (A, bufUniFactors, degree (A, y) + 1, Pi, diophant, M);
613    earlyFactors= earlyMonicFactorDetect (A, bufUniFactors, newLiftBound,
614                                           bufDegs, earlySuccess,
615                                           degree (A, y) + 1, liftBound);
616    if (bufDegs.getLength() > 1 && !earlySuccess)
617    {
618      if (newLiftBound > degree (A, y) + 1)
619      {
620        liftBound= newLiftBound;
621        bufUniFactors.insert (LC(A, x));
622        henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
623                            Pi, diophant, M);
624      }
625    }
626    else if (earlySuccess)
627      liftBound= newLiftBound;
628  }
629  else if (smallFactorDeg < degree (A, y) + 1)
630  {
631    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M);
632    earlyFactors= earlyMonicFactorDetect (A, bufUniFactors, newLiftBound,
633                                           bufDegs, earlySuccess,
634                                           smallFactorDeg, liftBound);
635    if (bufDegs.getLength() > 1 && !earlySuccess)
636    {
637      bufUniFactors.insert (LC (A, x));
638      henselLiftResume12 (A, bufUniFactors, smallFactorDeg, degree (A, y) +
639                          1, Pi, diophant, M);
640      earlyFactors= earlyMonicFactorDetect (A, bufUniFactors, newLiftBound,
641                                             bufDegs, earlySuccess,
642                                             degree (A, y) + 1, liftBound);
643      if (bufDegs.getLength() > 1 && !earlySuccess)
644      {
645        if (newLiftBound > degree (A, y) + 1)
646        {
647          if (newLiftBound < newLiftBound)
648            liftBound= newLiftBound;
649          bufUniFactors.insert (LC(A, x));
650          henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
651                              Pi, diophant, M);
652        }
653      }
654      else if (earlySuccess)
655        liftBound= newLiftBound;
656    }
657    else if (earlySuccess)
658      liftBound= newLiftBound;
659  }
660
661  if (newLiftBound > 0)
662    liftBound= newLiftBound;
663
664  CanonicalForm MODl= power (y, liftBound);
665
666  if (bufDegs.getLength() > 1)
667    factors= monicFactorRecombi (bufUniFactors, A, MODl, bufDegs);
668
669  if (earlySuccess)
670    factors= Union (earlyFactors, factors);
671  else if (!earlySuccess && bufDegs.getLength() == 1)
672    factors= earlyFactors;
673
674  decompressAppend (factors, contentAyFactors, N);
675
676  bivarEval= N (bivarEval);
677
678  return factors;
679}
680
681CFList
682extFactorRecombination (const CFList& factors, const CanonicalForm& F,
683                        const CFList& M, const ExtensionInfo& info,
684                        const CFList& evaluation)
685{
686  Variable alpha= info.getAlpha();
687  Variable beta= info.getBeta();
688  CanonicalForm gamma= info.getGamma();
689  CanonicalForm delta= info.getDelta();
690  int k= info.getGFDegree();
691  CFList source, dest;
692  if (factors.length() == 1)
693  {
694    CanonicalForm buf= reverseShift (F, evaluation);
695    return CFList (mapDown (buf, info, source, dest));
696  }
697  if (factors.length() < 1)
698    return CFList();
699
700  int degMipoBeta;
701  if (!k && beta == Variable(1))
702    degMipoBeta= 1;
703  else if (!k && beta != Variable(1))
704    degMipoBeta= degree (getMipo (beta));
705
706  CFList T, S;
707  T= factors;
708
709  int s= 1;
710  CFList result;
711  CanonicalForm buf;
712  if (beta != Variable (1))
713    buf= mapDown (F, gamma, delta, alpha, source, dest);
714  else
715    buf= F;
716
717  CanonicalForm g, LCBuf= LC (buf, Variable (1));
718  CanonicalForm buf2;
719  int v [T.length()];
720  for (int i= 0; i < T.length(); i++)
721    v[i]= 0;
722  bool noSubset= false;
723  CFArray TT;
724  int subsetDeg;
725  TT= copy (factors);
726  bool recombination= false;
727  while (T.length() >= 2*s)
728  {
729    while (noSubset == false)
730    {
731      if (T.length() == s)
732      {
733        if (recombination)
734        {
735          T.insert (LCBuf);
736          g= prodMod (T, M);
737          T.removeFirst();
738          result.append (g/myContent (g));
739          g= reverseShift (g, evaluation);
740          g /= Lc (g);
741          appendTestMapDown (result, g, info, source, dest);
742          return result;
743        }
744        else
745        {
746          buf= reverseShift (buf, evaluation);
747          return CFList (buf);
748        }
749        return result;
750      }
751
752      S= subset (v, s, TT, noSubset);
753      if (noSubset) break;
754
755      S.insert (LCBuf);
756      g= prodMod (S, M);
757      S.removeFirst();
758      g /= myContent (g);
759      if (fdivides (g, buf))
760      {
761        buf2= reverseShift (g, evaluation);
762        buf2 /= Lc (buf2);
763        if (!k && beta == Variable (1))
764        {
765          if (degree (buf2, alpha) < degMipoBeta)
766          {
767            appendTestMapDown (result, buf2, info, source, dest);
768            buf /= g;
769            LCBuf= LC (buf, Variable (1));
770            recombination= true;
771          }
772        }
773        else
774        {
775          if (!isInExtension (buf2, delta, k))
776          {
777            appendTestMapDown (result, buf2, info, source, dest);
778            buf /= g;
779            LCBuf= LC (buf, Variable (1));
780            recombination= true;
781          }
782        }
783        T= Difference (T, S);
784
785        if (T.length() < 2*s || T.length() == s)
786        {
787          buf= reverseShift (buf, evaluation);
788          buf /= Lc (buf);
789          appendTestMapDown (result, buf, info, source, dest);
790          return result;
791        }
792        TT= copy (T);
793        indexUpdate (v, s, T.length(), noSubset);
794        if (noSubset) break;
795      }
796    }
797    s++;
798    if (T.length() < 2*s || T.length() == s)
799    {
800      buf= reverseShift (buf, evaluation);
801      appendTestMapDown (result, buf, info, source, dest);
802      return result;
803    }
804    int v [T.length()];
805    for (int i= 0; i < T.length(); i++)
806      v[i]= 0;
807    noSubset= false;
808  }
809  if (T.length() < 2*s)
810  {
811    buf= reverseShift (F, evaluation);
812    appendMapDown (result, buf, info, source, dest);
813  }
814
815  return result;
816}
817
818CFList
819factorRecombination (const CanonicalForm& F, const CFList& factors,
820                     const CFList& M)
821{
822  if (factors.length() == 1)
823    return CFList(F);
824  if (factors.length() < 1)
825    return CFList();
826
827  CFList T, S;
828
829  T= factors;
830
831  int s= 1;
832  CFList result;
833  CanonicalForm LCBuf= LC (F, Variable (1));
834  CanonicalForm g, buf= F;
835  int v [T.length()];
836  for (int i= 0; i < T.length(); i++)
837    v[i]= 0;
838  bool noSubset= false;
839  CFArray TT;
840  int subsetDeg;
841  TT= copy (factors);
842  Variable y= F.level() - 1;
843  bool recombination= false;
844  CanonicalForm h;
845  while (T.length() >= 2*s)
846  {
847    while (noSubset == false)
848    {
849      if (T.length() == s)
850      {
851        if (recombination)
852        {
853          T.insert (LC (buf));
854          g= prodMod (T, M);
855          result.append (g/myContent (g));
856          return result;
857        }
858        else
859          return CFList (F);
860      }
861      S= subset (v, s, TT, noSubset);
862      if (noSubset) break;
863      S.insert (LCBuf);
864      g= prodMod (S, M);
865      S.removeFirst();
866      g /= myContent (g);
867      if (fdivides (g, buf))
868      {
869        recombination= true;
870        result.append (g);
871        buf /= g;
872        LCBuf= LC (buf, Variable(1));
873        T= Difference (T, S);
874        if (T.length() < 2*s || T.length() == s)
875        {
876          result.append (buf);
877          return result;
878        }
879        TT= copy (T);
880        indexUpdate (v, s, T.length(), noSubset);
881        if (noSubset) break;
882      }
883    }
884    s++;
885    if (T.length() < 2*s || T.length() == s)
886    {
887      result.append (buf);
888      return result;
889    }
890    int v [T.length()];
891    for (int i= 0; i < T.length(); i++)
892      v[i]= 0;
893    noSubset= false;
894  }
895  if (T.length() < 2*s)
896    result.append (F);
897
898  return result;
899}
900
901int
902liftBoundAdaption (const CanonicalForm& F, const CFList& factors, bool&
903                   success, const int deg, const CFList& MOD, const int bound)
904{
905  int adaptedLiftBound= 0;
906  CanonicalForm buf= F;
907  Variable y= F.mvar();
908  CanonicalForm LCBuf= LC (buf, Variable (1));
909  CanonicalForm g;
910  CFList M= MOD;
911  M.append (power (y, deg));
912  int d= bound;
913  int e= 0;
914  int nBuf;
915  for (CFListIterator i= factors; i.hasItem(); i++)
916  {
917    g= mulMod (i.getItem(), LCBuf, M);
918    g /= myContent (g);
919    if (fdivides (g, buf))
920    {
921      nBuf= degree (g, y) + degree (LC (g, 1), y);
922      d -= nBuf;
923      e= tmax (e, nBuf);
924      buf /= g;
925      LCBuf= LC (buf, Variable (1));
926    }
927  }
928  adaptedLiftBound= d;
929
930  if (adaptedLiftBound < deg)
931  {
932    if (adaptedLiftBound < degree (F) + 1)
933    {
934      if (d == 1)
935      {
936        if (e + 1 > deg)
937        {
938          adaptedLiftBound= deg;
939          success= false;
940        }
941        else
942        {
943          success= true;
944          if (e + 1 < degree (F) + 1)
945            adaptedLiftBound= deg;
946          else
947            adaptedLiftBound= e + 1;
948        }
949      }
950      else
951      {
952        success= true;
953        adaptedLiftBound= deg;
954      }
955    }
956    else
957    {
958      success= true;
959    }
960  }
961  return adaptedLiftBound;
962}
963
964int
965extLiftBoundAdaption (const CanonicalForm& F, const CFList& factors, bool&
966                      success, const ExtensionInfo& info, const CFList& eval,
967                      const int deg, const CFList& MOD, const int bound)
968{
969  Variable alpha= info.getAlpha();
970  Variable beta= info.getBeta();
971  CanonicalForm gamma= info.getGamma();
972  CanonicalForm delta= info.getDelta();
973  int k= info.getGFDegree();
974  int adaptedLiftBound= 0;
975  CanonicalForm buf= F;
976  Variable y= F.mvar();
977  CanonicalForm LCBuf= LC (buf, Variable (1));
978  CanonicalForm g, gg;
979  CFList M= MOD;
980  M.append (power (y, deg));
981  adaptedLiftBound= 0;
982  int d= bound;
983  int e= 0;
984  int nBuf;
985  int degMipoBeta;
986  if (!k && beta == Variable(1))
987    degMipoBeta= 1;
988  else if (!k && beta != Variable(1))
989    degMipoBeta= degree (getMipo (beta));
990
991  for (CFListIterator i= factors; i.hasItem(); i++)
992  {
993    g= mulMod (i.getItem(), LCBuf, M);
994    g /= myContent (g);
995    if (fdivides (g, buf))
996    {
997      gg= reverseShift (g, eval);
998      gg /= Lc (gg);
999      if (!k && beta == Variable (1))
1000      {
1001        if (degree (gg, alpha) < degMipoBeta)
1002        {
1003          buf /= g;
1004          nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1005          d -= nBuf;
1006          e= tmax (e, nBuf);
1007          LCBuf= LC (buf, Variable (1));
1008        }
1009      }
1010      else
1011      {
1012        if (!isInExtension (gg, delta, k))
1013        {
1014          buf /= g;
1015          nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1016          d -= nBuf;
1017          e= tmax (e, nBuf);
1018          LCBuf= LC (buf, Variable (1));
1019        }
1020      }
1021    }
1022  }
1023  adaptedLiftBound= d;
1024
1025  if (adaptedLiftBound < deg)
1026  {
1027    if (adaptedLiftBound < degree (F) + 1)
1028    {
1029      if (d == 1)
1030      {
1031        if (e + 1 > deg)
1032        {
1033          adaptedLiftBound= deg;
1034          success= false;
1035        }
1036        else
1037        {
1038          success= true;
1039          if (e + 1 < degree (F) + 1)
1040            adaptedLiftBound= deg;
1041          else
1042            adaptedLiftBound= e + 1;
1043        }
1044      }
1045      else
1046      {
1047        success= true;
1048        adaptedLiftBound= deg;
1049      }
1050    }
1051    else
1052    {
1053      success= true;
1054    }
1055  }
1056
1057  return adaptedLiftBound;
1058}
1059
1060CFList
1061earlyFactorDetect (CanonicalForm& F, CFList& factors, int& adaptedLiftBound,
1062                   bool& success, const int deg, const CFList& MOD,
1063                   const int bound)
1064{
1065  CFList result;
1066  CFList T= factors;
1067  CanonicalForm buf= F;
1068  Variable y= F.mvar();
1069  CanonicalForm LCBuf= LC (buf, Variable (1));
1070  CanonicalForm g;
1071  CFList M= MOD;
1072  M.append (power (y, deg));
1073  adaptedLiftBound= 0;
1074  int d= bound;
1075  int e= 0;
1076  int nBuf;
1077  for (CFListIterator i= factors; i.hasItem(); i++)
1078  {
1079    g= mulMod (i.getItem(), LCBuf, M);
1080    g /= myContent (g);
1081    if (fdivides (g, buf))
1082    {
1083      result.append (g);
1084      nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1085      d -= nBuf;
1086      e= tmax (e, nBuf);
1087      buf /= g;
1088      LCBuf= LC (buf, Variable (1));
1089      T= Difference (T, CFList (i.getItem()));
1090    }
1091  }
1092  adaptedLiftBound= d;
1093
1094  if (adaptedLiftBound < deg)
1095  {
1096    if (adaptedLiftBound < degree (F) + 1)
1097    {
1098      if (d == 1)
1099        adaptedLiftBound= tmin (e + 1, deg);
1100      else
1101        adaptedLiftBound= deg;
1102    }
1103    factors= T;
1104    F= buf;
1105    success= true;
1106  }
1107  return result;
1108}
1109
1110CFList
1111extEarlyFactorDetect (CanonicalForm& F, CFList& factors, int& adaptedLiftBound,
1112                      bool& success, const ExtensionInfo& info, const CFList&
1113                      eval, const int deg, const CFList& MOD, const int bound)
1114{
1115  Variable alpha= info.getAlpha();
1116  Variable beta= info.getBeta();
1117  CanonicalForm gamma= info.getGamma();
1118  CanonicalForm delta= info.getDelta();
1119  int k= info.getGFDegree();
1120  CFList result;
1121  CFList T= factors;
1122  CanonicalForm buf= F;
1123  Variable y= F.mvar();
1124  CanonicalForm LCBuf= LC (buf, Variable (1));
1125  CanonicalForm g, gg;
1126  CFList M= MOD;
1127  M.append (power (y, deg));
1128  adaptedLiftBound= 0;
1129  int d= bound;
1130  int e= 0;
1131  int nBuf;
1132  CFList source, dest;
1133
1134  int degMipoBeta;
1135  if (!k && beta == Variable(1))
1136    degMipoBeta= 1;
1137  else if (!k && beta != Variable(1))
1138    degMipoBeta= degree (getMipo (beta));
1139
1140  for (CFListIterator i= factors; i.hasItem(); i++)
1141  {
1142    g= mulMod (i.getItem(), LCBuf, M);
1143    g /= myContent (g);
1144    if (fdivides (g, buf))
1145    {
1146      gg= reverseShift (g, eval);
1147      gg /= Lc (gg);
1148          if (!k && beta == Variable (1))
1149          {
1150            if (degree (gg, alpha) < degMipoBeta)
1151            {
1152              appendTestMapDown (result, gg, info, source, dest);
1153              buf /= g;
1154              nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1155              d -= nBuf;
1156              e= tmax (e, nBuf);
1157              LCBuf= LC (buf, Variable (1));
1158            }
1159          }
1160          else
1161          {
1162            if (!isInExtension (gg, delta, k))
1163            {
1164              appendTestMapDown (result, gg, info, source, dest);
1165              buf /= g;
1166              nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1167              d -= nBuf;
1168              e= tmax (e, nBuf);
1169              LCBuf= LC (buf, Variable (1));
1170            }
1171          }
1172      T= Difference (T, CFList (i.getItem()));
1173    }
1174  }
1175  adaptedLiftBound= d;
1176
1177  if (adaptedLiftBound < deg)
1178  {
1179    if (adaptedLiftBound < degree (F) + 1)
1180    {
1181      if (d == 1)
1182        adaptedLiftBound= tmin (e + 1, deg);
1183      else
1184        adaptedLiftBound= deg;
1185    }
1186    success= true;
1187    factors= T;
1188    F= buf;
1189  }
1190  return result;
1191}
1192
1193CFList
1194evalPoints (const CanonicalForm& F, CFList & eval, const Variable& alpha,
1195            CFList& list, const bool& GF, bool& fail)
1196{
1197  int k= F.level() - 1;
1198  Variable x= Variable (1);
1199  CFList result;
1200  FFRandom genFF;
1201  GFRandom genGF;
1202  int p= getCharacteristic ();
1203  double bound;
1204  if (alpha != Variable (1))
1205  {
1206    bound= pow ((double) p, (double) degree (getMipo(alpha)));
1207    bound= pow ((double) bound, (double) k);
1208  }
1209  else if (GF)
1210  {
1211    bound= pow ((double) p, (double) getGFDegree());
1212    bound= pow ((double) bound, (double) k);
1213  }
1214  else
1215    bound= pow ((double) p, (double) k);
1216
1217  CanonicalForm random;
1218  CanonicalForm deriv_x, gcd_deriv;
1219  do
1220  {
1221    random= 0;
1222    // possible overflow if list.length() does not fit into a int
1223    if (list.length() >= bound)
1224    {
1225      fail= true;
1226      break;
1227    }
1228    for (int i= 0; i < k; i++)
1229    {
1230      if (GF)
1231      {
1232        result.append (genGF.generate());
1233        random += result.getLast()*power (x, i);
1234      }
1235      else if (alpha != Variable(1))
1236      {
1237        AlgExtRandomF genAlgExt (alpha);
1238        result.append (genAlgExt.generate());
1239        random += result.getLast()*power (x, i);
1240      }
1241      else
1242      {
1243        result.append (genFF.generate());
1244        random += result.getLast()*power (x, i);
1245      }
1246    }
1247    if (find (list, random))
1248    {
1249      result= CFList();
1250      continue;
1251    }
1252    int l= F.level();
1253    eval.insert (F);
1254    for (CFListIterator i= result; i.hasItem(); i++, l--)
1255      eval.insert (eval.getFirst()(i.getItem(), l));
1256
1257    if (degree (eval.getFirst()) != degree (F, 1))
1258    {
1259      if (!find (list, random))
1260        list.append (random);
1261      result= CFList();
1262      eval= CFList();
1263      continue;
1264    }
1265
1266    deriv_x= deriv (eval.getFirst(), x);
1267    gcd_deriv= gcd (eval.getFirst(), deriv_x);
1268    if (degree (gcd_deriv) > 0)
1269    {
1270      if (!find (list, random))
1271        list.append (random);
1272      result= CFList();
1273      eval= CFList();
1274      continue;
1275    }
1276    CFListIterator i= eval;
1277    i++;
1278    CanonicalForm contentx= content (i.getItem(), x);
1279    if (degree (contentx) > 0)
1280    {
1281      if (!find (list, random))
1282        list.append (random);
1283      result= CFList();
1284      eval= CFList();
1285      continue;
1286    }
1287
1288    if (list.length() >= bound)
1289    {
1290      fail= true;
1291      break;
1292    }
1293  } while (find (list, random));
1294
1295  if (!eval.isEmpty())
1296    eval.removeFirst();
1297
1298  return result;
1299}
1300
1301static inline
1302int newMainVariableSearch (CanonicalForm& A, CFList& Aeval, CFList&
1303                           evaluation, const Variable& alpha, const int lev)
1304{
1305  Variable x= Variable (1);
1306  CanonicalForm derivI, buf;
1307  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1308  int swapLevel= 0;
1309  CFList list;
1310  bool fail= false;
1311  buf= A;
1312  Aeval= CFList();
1313  evaluation= CFList();
1314  for (int i= lev; i <= A.level(); i++)
1315  {
1316    derivI= deriv (buf, Variable (i));
1317    if (!derivI.isZero())
1318    {
1319      buf= swapvar (buf, x, Variable (i));
1320      Aeval= CFList();
1321      evaluation= CFList();
1322      fail= false;
1323      evaluation= evalPoints (buf, Aeval, alpha, list, GF, fail);
1324      if (!fail)
1325      {
1326        A= buf;
1327        swapLevel= i;
1328        break;
1329      }
1330      else
1331        buf= A;
1332    }
1333  }
1334  return swapLevel;
1335}
1336
1337static inline
1338CanonicalForm lcmContent (const CanonicalForm& A, CFList& contentAi)
1339{
1340  int i= A.level();
1341  contentAi.append (myContent (A, i));
1342  contentAi.append (myContent (A, i - 1));
1343  CanonicalForm result= myLcm (contentAi.getFirst(), contentAi.getLast());
1344  for (i= i - 2; i > 0; i--)
1345  {
1346    contentAi.append (content (A, i));
1347    result= myLcm (result, contentAi.getLast());
1348  }
1349  return result;
1350}
1351
1352CFList
1353henselLiftAndEarly (CanonicalForm& A, CFList& MOD, int*& liftBounds, bool&
1354                    earlySuccess, CFList& earlyFactors, const CFList& Aeval,
1355                    const CFList& biFactors, const CFList& evaluation,
1356                    const ExtensionInfo& info)
1357{
1358  bool extension= info.isInExtension();
1359  CFList bufFactors= biFactors;
1360  bufFactors.insert (LC (Aeval.getFirst(), 1));
1361
1362  sortList (bufFactors, Variable (1));
1363
1364  CFList diophant;
1365  CFArray Pi;
1366  int smallFactorDeg= 11; //tunable parameter
1367  CFList result;
1368  int newLiftBound= 0;
1369  int adaptedLiftBound= 0;
1370  int liftBound= liftBounds[1];
1371
1372  earlySuccess= false;
1373  bool earlyReconst= false;
1374  CFList earlyReconstFactors;
1375  CFListIterator j= Aeval;
1376  j++;
1377  CanonicalForm buf= j.getItem();
1378  CFMatrix Mat= CFMatrix (liftBound, bufFactors.length() - 1);
1379  MOD= CFList (power (Variable (2), liftBounds[0]));
1380  if (smallFactorDeg >= liftBound)
1381  {
1382    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
1383  }
1384  else if (smallFactorDeg >= degree (buf) + 1)
1385  {
1386    liftBounds[1]= degree (buf) + 1;
1387    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
1388    if (Aeval.length() == 2)
1389    {
1390      if (!extension)
1391        earlyFactors= earlyFactorDetect
1392                       (buf, result, adaptedLiftBound, earlySuccess,
1393                        degree (buf) + 1, MOD, liftBound);
1394      else
1395        earlyFactors= extEarlyFactorDetect
1396                       (buf, result, adaptedLiftBound, earlySuccess,
1397                        info, evaluation, degree
1398                        (buf) + 1, MOD, liftBound);
1399    }
1400    else
1401    {
1402      if (!extension)
1403        adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
1404                                             degree (buf) + 1, MOD, liftBound);
1405      else
1406        adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess, info,
1407                                                evaluation, degree (buf) + 1,
1408                                                MOD, liftBound);
1409    }
1410    if (!earlySuccess)
1411    {
1412      result.insert (LC (buf, 1));
1413      liftBounds[1]= adaptedLiftBound;
1414      liftBound= adaptedLiftBound;
1415      henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1416                        Pi, diophant, Mat, MOD);
1417    }
1418    else
1419      liftBounds[1]= adaptedLiftBound;
1420  }
1421  else if (smallFactorDeg < degree (buf) + 1)
1422  {
1423    liftBounds[1]= smallFactorDeg;
1424    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
1425    if (Aeval.length() == 2)
1426    {
1427      if (!extension)
1428        earlyFactors= earlyFactorDetect (buf, result, adaptedLiftBound,
1429                                         earlySuccess, smallFactorDeg, MOD,
1430                                         liftBound);
1431      else
1432        earlyFactors= extEarlyFactorDetect (buf, result, adaptedLiftBound,
1433                                            earlySuccess, info, evaluation,
1434                                            smallFactorDeg, MOD, liftBound);
1435    }
1436    else
1437    {
1438      if (!extension)
1439        adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
1440                                             smallFactorDeg, MOD, liftBound);
1441      else
1442        adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess, info,
1443                                                evaluation, smallFactorDeg, MOD,
1444                                                liftBound);
1445    }
1446
1447    if (!earlySuccess)
1448    {
1449      result.insert (LC (buf, 1));
1450      henselLiftResume (buf, result, smallFactorDeg, degree (buf) + 1,
1451                        Pi, diophant, Mat, MOD);
1452      if (Aeval.length() == 2)
1453      {
1454         if (!extension)
1455           earlyFactors= earlyFactorDetect (buf, result, adaptedLiftBound,
1456                                            earlySuccess, degree (buf) + 1,
1457                                            MOD, liftBound);
1458         else
1459           earlyFactors= extEarlyFactorDetect (buf, result, adaptedLiftBound,
1460                                               earlySuccess, info, evaluation,
1461                                               degree (buf) + 1, MOD,
1462                                               liftBound);
1463      }
1464      else
1465      {
1466        if (!extension)
1467          adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
1468                                               degree (buf) + 1, MOD,liftBound);
1469        else
1470          adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess,
1471                                                  info, evaluation,
1472                                                  degree (buf) + 1, MOD,
1473                                                  liftBound);
1474      }
1475      if (!earlySuccess)
1476      {
1477        result.insert (LC (buf, 1));
1478        liftBounds[1]= adaptedLiftBound;
1479        liftBound= adaptedLiftBound;
1480        henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1481                          Pi, diophant, Mat, MOD);
1482      }
1483      else
1484        liftBounds[1]= adaptedLiftBound;
1485    }
1486    else
1487      liftBounds[1]= adaptedLiftBound;
1488  }
1489
1490  MOD.append (power (Variable (3), liftBounds[1]));
1491
1492  if (Aeval.length() > 2)
1493  {
1494    CFListIterator j= Aeval;
1495    j++;
1496    CFList bufEval;
1497    bufEval.append (j.getItem());
1498    j++;
1499    int liftBoundsLength= Aeval.getLast().level() - 1;
1500    for (int i= 2; i <= liftBoundsLength && j.hasItem(); i++, j++)
1501    {
1502      earlySuccess= false;
1503      result.insert (LC (bufEval.getFirst(), 1));
1504      bufEval.append (j.getItem());
1505      liftBound= liftBounds[i];
1506      Mat= CFMatrix (liftBounds[i], result.length() - 1);
1507
1508      buf= j.getItem();
1509      if (smallFactorDeg >= liftBound)
1510        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
1511                            liftBounds[i -  1], liftBounds[i]);
1512      else if (smallFactorDeg >= degree (buf) + 1)
1513      {
1514        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
1515                            liftBounds[i -  1], degree (buf) + 1);
1516
1517        if (Aeval.length() == i + 1)
1518        {
1519          if (!extension)
1520            earlyFactors= earlyFactorDetect
1521                           (buf, result, adaptedLiftBound, earlySuccess,
1522                            degree (buf) + 1, MOD, liftBound);
1523          else
1524            earlyFactors= extEarlyFactorDetect
1525                           (buf, result, adaptedLiftBound, earlySuccess,
1526                            info, evaluation, degree (buf) + 1, MOD, liftBound);
1527        }
1528        else
1529        {
1530          if (!extension)
1531            adaptedLiftBound= liftBoundAdaption
1532                                (buf, result, earlySuccess, degree (buf)
1533                                 + 1,  MOD, liftBound);
1534          else
1535            adaptedLiftBound= extLiftBoundAdaption
1536                                (buf, result, earlySuccess, info, evaluation,
1537                                 degree (buf) + 1, MOD, liftBound);
1538        }
1539
1540        if (!earlySuccess)
1541        {
1542          result.insert (LC (buf, 1));
1543          liftBounds[i]= adaptedLiftBound;
1544          liftBound= adaptedLiftBound;
1545          henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1546                            Pi, diophant, Mat, MOD);
1547        }
1548        else
1549        {
1550          liftBounds[i]= adaptedLiftBound;
1551        }
1552      }
1553      else if (smallFactorDeg < degree (buf) + 1)
1554      {
1555        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
1556                            liftBounds[i -  1], smallFactorDeg);
1557
1558        if (Aeval.length() == i + 1)
1559        {
1560          if (!extension)
1561            earlyFactors= earlyFactorDetect
1562                           (buf, result, adaptedLiftBound, earlySuccess,
1563                            smallFactorDeg, MOD, liftBound);
1564          else
1565            earlyFactors= extEarlyFactorDetect
1566                           (buf, result, adaptedLiftBound, earlySuccess,
1567                            info, evaluation, smallFactorDeg, MOD, liftBound);
1568        }
1569        else
1570        {
1571          if (!extension)
1572            adaptedLiftBound= liftBoundAdaption
1573                                (buf, result, earlySuccess,
1574                                 smallFactorDeg, MOD, liftBound);
1575          else
1576            adaptedLiftBound= extLiftBoundAdaption
1577                                (buf, result, earlySuccess, info, evaluation,
1578                                 smallFactorDeg, MOD, liftBound);
1579        }
1580
1581        if (!earlySuccess)
1582        {
1583          result.insert (LC (buf, 1));
1584          henselLiftResume (buf, result, smallFactorDeg,
1585                            degree (buf) + 1, Pi, diophant, Mat, MOD);
1586          if (Aeval.length() == i + 1)
1587          {
1588            if (!extension)
1589              earlyFactors= earlyFactorDetect
1590                             (buf, result, adaptedLiftBound, earlySuccess,
1591                              degree (buf) +  1,  MOD, liftBound);
1592            else
1593              earlyFactors= extEarlyFactorDetect
1594                             (buf, result, adaptedLiftBound, earlySuccess,
1595                              info, evaluation, degree (buf) + 1, MOD,
1596                              liftBound);
1597          }
1598          else
1599          {
1600            if (!extension)
1601              adaptedLiftBound= liftBoundAdaption
1602                                  (buf, result, earlySuccess, degree
1603                                   (buf) +  1,  MOD, liftBound);
1604            else
1605              adaptedLiftBound= extLiftBoundAdaption
1606                                  (buf, result, earlySuccess, info, evaluation,
1607                                   degree (buf) + 1,  MOD, liftBound);
1608          }
1609
1610          if (!earlySuccess)
1611          {
1612            result.insert (LC (buf, 1));
1613            liftBounds[i]= adaptedLiftBound;
1614            liftBound= adaptedLiftBound;
1615            henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1616                              Pi, diophant, Mat, MOD);
1617          }
1618          else
1619            liftBounds[i]= adaptedLiftBound;
1620        }
1621        else
1622          liftBounds[i]= adaptedLiftBound;
1623      }
1624      MOD.append (power (Variable (i + 2), liftBounds[i]));
1625      bufEval.removeFirst();
1626    }
1627    bufFactors= result;
1628  }
1629  else
1630    bufFactors= result;
1631
1632  if (earlySuccess)
1633    A= buf;
1634  return result;
1635}
1636
1637CFList
1638extFactorize (const CanonicalForm& F, const ExtensionInfo& info);
1639
1640CFList
1641multiFactorize (const CanonicalForm& F, const ExtensionInfo& info)
1642{
1643
1644  if (F.inCoeffDomain())
1645    return CFList (F);
1646
1647  // compress and find main Variable
1648  CFMap N;
1649  CanonicalForm A= myCompress (F, N);
1650
1651  A /= Lc (A); // make monic
1652
1653  Variable alpha= info.getAlpha();
1654  Variable beta= info.getBeta();
1655  CanonicalForm gamma= info.getGamma();
1656  CanonicalForm delta= info.getDelta();
1657  int k= info.getGFDegree();
1658  bool extension= info.isInExtension();
1659  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1660  //univariate case
1661  if (F.isUnivariate())
1662  {
1663    if (extension == false)
1664      return uniFactorizer (F, alpha, GF);
1665    else
1666    {
1667      CFList source, dest;
1668      A= mapDown (F, info, source, dest);
1669      return uniFactorizer (A, beta, GF);
1670    }
1671  }
1672
1673  //bivariate case
1674  if (A.level() == 2)
1675  {
1676    CFList buf= biFactorize (F, info);
1677    return buf;
1678  }
1679
1680  Variable x= Variable (1);
1681  Variable y= Variable (2);
1682
1683  // remove content
1684  CFList contentAi;
1685  CanonicalForm lcmCont= lcmContent (A, contentAi);
1686  A /= lcmCont;
1687
1688  // trivial after content removal
1689  CFList contentAFactors;
1690  if (A.inCoeffDomain())
1691  {
1692    for (CFListIterator i= contentAi; i.hasItem(); i++)
1693    {
1694      if (i.getItem().inCoeffDomain())
1695        continue;
1696      else
1697      {
1698        lcmCont /= i.getItem();
1699        contentAFactors=
1700        Union (multiFactorize (lcmCont, info),
1701               multiFactorize (i.getItem(), info));
1702        break;
1703      }
1704    }
1705    decompress (contentAFactors, N);
1706    normalize (contentAFactors);
1707    return contentAFactors;
1708  }
1709
1710  // factorize content
1711  contentAFactors= multiFactorize (lcmCont, info);
1712
1713  // univariate after content removal
1714  CFList factors;
1715  if (A.isUnivariate ())
1716  {
1717    factors= uniFactorizer (A, alpha, GF);
1718    append (factors, contentAFactors);
1719    decompress (factors, N);
1720    return factors;
1721  }
1722
1723  // check main variable
1724  int swapLevel= 0;
1725  CanonicalForm derivZ;
1726  CanonicalForm gcdDerivZ;
1727  CanonicalForm bufA= A;
1728  Variable z;
1729  for (int i= 1; i <= A.level(); i++)
1730  {
1731    z= Variable (i);
1732    derivZ= deriv (bufA, z);
1733    if (derivZ.isZero())
1734    {
1735      if (i == 1)
1736        swapLevel= 1;
1737      else
1738        continue;
1739    }
1740    else
1741    {
1742      if (swapLevel == 1)
1743      {
1744        swapLevel= i;
1745        A= swapvar (A, x, z);
1746      }
1747      if (GF)
1748        gcdDerivZ= GCD_GF (bufA, derivZ);
1749      else if (alpha == Variable (1))
1750        gcdDerivZ= GCD_small_p (bufA, derivZ);
1751      else
1752        gcdDerivZ= GCD_Fp_extension (bufA, derivZ, alpha);
1753      if (degree (gcdDerivZ) > 0 && !derivZ.isZero())
1754      {
1755        CanonicalForm g= bufA/gcdDerivZ;
1756        CFList factorsG=
1757        Union (multiFactorize (g, info),
1758               multiFactorize (gcdDerivZ, info));
1759        appendSwapDecompress (factorsG, contentAFactors, N, swapLevel, x);
1760        normalize (factorsG);
1761        return factorsG;
1762      }
1763    }
1764  }
1765
1766
1767  CFList Aeval, list, evaluation, bufEvaluation, bufAeval;
1768  bool fail= false;
1769  int swapLevel2= 0;
1770  int level;
1771  int factorNums= 3;
1772  CanonicalForm bivarEval;
1773  CFList biFactors, bufBiFactors;
1774  CanonicalForm evalPoly;
1775  int lift, bufLift;
1776  if (factorNums < (int) ilog2 (totaldegree (A)))
1777    factorNums= (int) ilog2 (totaldegree (A));
1778  // several bivariate factorizations
1779  for (int i= 0; i < factorNums; i++)
1780  {
1781    bufA= A;
1782    bufAeval= CFList();
1783    bufEvaluation= evalPoints (bufA, bufAeval, alpha, list, GF, fail);
1784    evalPoly= 0;
1785
1786    if (fail && (i == 0))
1787    {
1788      if (!swapLevel)
1789        level= 2;
1790      else
1791        level= swapLevel + 1;
1792
1793      swapLevel2= newMainVariableSearch (A, Aeval, evaluation, alpha, level);
1794
1795      if (!swapLevel2) // need to pass to an extension
1796      {
1797        factors= extFactorize (A, info);
1798        appendSwapDecompress (factors, contentAFactors, N, swapLevel, x);
1799        normalize (factors);
1800        return factors;
1801      }
1802      else
1803      {
1804        fail= false;
1805        bufAeval= Aeval;
1806        bufA= A;
1807        bufEvaluation= evaluation;
1808      }
1809    }
1810    else if (fail && (i > 0))
1811      break;
1812
1813    bivarEval= bufEvaluation.getLast();
1814    bufEvaluation.removeLast();
1815
1816    bufLift= degree (A, y) + 1 + degree (LC(A, x), y);
1817
1818    TIMING_START (fac_bi_factorizer);
1819    bufBiFactors= biFactorizer (bufAeval.getFirst(), alpha, bivarEval, bufLift);
1820    TIMING_END_AND_PRINT (fac_bi_factorizer,
1821                          "time for bivariate factorization: ");
1822
1823    if (bufBiFactors.length() == 1)
1824    {
1825      if (extension)
1826      {
1827        CFList source, dest;
1828        A= mapDown (A, info, source, dest);
1829      }
1830      factors.append (A);
1831      appendSwapDecompress (factors, contentAFactors, N, swapLevel,
1832                            swapLevel2, x);
1833      normalize (factors);
1834      return factors;
1835    }
1836
1837    bufEvaluation.append (-bivarEval[0]);
1838    if (i == 0)
1839    {
1840      Aeval= bufAeval;
1841      evaluation= bufEvaluation;
1842      biFactors= bufBiFactors;
1843      lift= bufLift;
1844    }
1845    else
1846    {
1847      if (bufBiFactors.length() < biFactors.length() ||
1848          ((bufLift < lift) && (bufBiFactors.length() == biFactors.length())))
1849      {
1850        Aeval= bufAeval;
1851        evaluation= bufEvaluation;
1852        biFactors= bufBiFactors;
1853        lift= bufLift;
1854      }
1855    }
1856    int k= 0;
1857    for (CFListIterator j= bufEvaluation; j.hasItem(); j++, k++)
1858      evalPoly += j.getItem()*power (x, k);
1859    list.append (evalPoly);
1860  }
1861
1862  //shifting to zero
1863  A= shift2Zero (A, Aeval, evaluation);
1864
1865  int* liftBounds;
1866  int liftBoundsLength= F.level() - 1;
1867  liftBounds= liftingBounds (A, lift);
1868
1869  CFList MOD;
1870  bool earlySuccess;
1871  CFList earlyFactors;
1872  TIMING_START (fac_hensel_lift);
1873  CFList liftedFactors= henselLiftAndEarly
1874                        (A, MOD, liftBounds, earlySuccess, earlyFactors,
1875                         Aeval, biFactors, evaluation, info);
1876  TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
1877
1878  if (!extension)
1879  {
1880    TIMING_START (fac_factor_recombination);
1881    factors= factorRecombination (A, liftedFactors, MOD);
1882    TIMING_END_AND_PRINT (fac_factor_recombination,
1883                          "time for factor recombination: ");
1884  }
1885  else
1886  {
1887    TIMING_START (fac_factor_recombination);
1888    factors= extFactorRecombination (liftedFactors, A, MOD, info, evaluation);
1889    TIMING_END_AND_PRINT (fac_factor_recombination,
1890                          "time for factor recombination: ");
1891  }
1892
1893  if (earlySuccess)
1894    factors= Union (factors, earlyFactors);
1895
1896  if (!extension)
1897  {
1898    for (CFListIterator i= factors; i.hasItem(); i++)
1899    {
1900      int kk= Aeval.getLast().level();
1901      for (CFListIterator j= evaluation; j.hasItem(); j++, kk--)
1902      {
1903        if (i.getItem().level() < kk)
1904          continue;
1905        i.getItem()= i.getItem() (Variable (kk) - j.getItem(), kk);
1906      }
1907    }
1908  }
1909
1910  swap (factors, swapLevel, swapLevel2, x);
1911  append (factors, contentAFactors);
1912  decompress (factors, N);
1913  normalize (factors);
1914
1915  delete[] liftBounds;
1916
1917  return factors;
1918}
1919
1920/// multivariate factorization over an extension of the initial field
1921CFList
1922extFactorize (const CanonicalForm& F, const ExtensionInfo& info)
1923{
1924  CanonicalForm A= F;
1925
1926  Variable alpha= info.getAlpha();
1927  Variable beta= info.getBeta();
1928  int k= info.getGFDegree();
1929  char cGFName= info.getGFName();
1930  CanonicalForm delta= info.getDelta();
1931  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1932  Variable w= Variable (1);
1933
1934  CFList factors;
1935  if (!GF && alpha == w)  // we are in F_p
1936  {
1937    CFList factors;
1938    bool extension= true;
1939    int p= getCharacteristic();
1940    if (p*p < (1<<16)) // pass to GF if possible
1941    {
1942      setCharacteristic (getCharacteristic(), 2, 'Z');
1943      ExtensionInfo info= ExtensionInfo (extension);
1944      A= A.mapinto();
1945      factors= multiFactorize (A, info);
1946
1947      Variable vBuf= rootOf (gf_mipo);
1948      setCharacteristic (getCharacteristic());
1949      for (CFListIterator j= factors; j.hasItem(); j++)
1950        j.getItem()= GF2FalphaRep (j.getItem(), vBuf);
1951    }
1952    else  // not able to pass to GF, pass to F_p(\alpha)
1953    {
1954      Variable v= chooseExtension (A, beta);
1955      ExtensionInfo info= ExtensionInfo (v, extension);
1956      factors= multiFactorize (A, info);
1957    }
1958    return factors;
1959  }
1960  else if (!GF && (alpha != w)) // we are in F_p(\alpha)
1961  {
1962    if (k == 1) // need factorization over F_p
1963    {
1964      int extDeg= degree (getMipo (alpha));
1965      extDeg++;
1966      CanonicalForm mipo= randomIrredpoly (extDeg + 1, Variable (1));
1967      Variable v= rootOf (mipo);
1968      ExtensionInfo info= ExtensionInfo (v);
1969      factors= biFactorize (A, info);
1970    }
1971    else
1972    {
1973      if (beta == Variable (1))
1974      {
1975        Variable v= chooseExtension (A, alpha);
1976        CanonicalForm primElem, imPrimElem;
1977        bool primFail= false;
1978        Variable vBuf;
1979        primElem= primitiveElement (alpha, vBuf, primFail);
1980        ASSERT (!primFail, "failure in integer factorizer");
1981        if (primFail)
1982          ; //ERROR
1983        else
1984          imPrimElem= mapPrimElem (primElem, vBuf, v);
1985
1986        CFList source, dest;
1987        CanonicalForm bufA= mapUp (A, alpha, v, primElem, imPrimElem,
1988                                   source, dest);
1989        ExtensionInfo info= ExtensionInfo (v, alpha, imPrimElem, primElem);
1990        factors= biFactorize (bufA, info);
1991      }
1992      else
1993      {
1994        Variable v= chooseExtension (A, alpha);
1995        CanonicalForm primElem, imPrimElem;
1996        bool primFail= false;
1997        Variable vBuf;
1998        ASSERT (!primFail, "failure in integer factorizer");
1999        if (primFail)
2000          ; //ERROR
2001        else
2002          imPrimElem= mapPrimElem (delta, beta, v); //oder mapPrimElem (primElem, vBuf, v);
2003
2004        CFList source, dest;
2005        CanonicalForm bufA= mapDown (A, info, source, dest);
2006        source= CFList();
2007        dest= CFList();
2008        bufA= mapUp (bufA, beta, v, delta, imPrimElem, source, dest);
2009        ExtensionInfo info= ExtensionInfo (v, beta, imPrimElem, delta);
2010        factors= biFactorize (bufA, info);
2011      }
2012    }
2013    return factors;
2014  }
2015  else // we are in GF (p^k)
2016  {
2017    int p= getCharacteristic();
2018    int extensionDeg= getGFDegree();
2019    bool extension= true;
2020    if (k == 1) // need factorization over F_p
2021    {
2022      extensionDeg++;
2023      if (pow ((double) p, (double) extensionDeg) < (1<<16))
2024      // pass to GF(p^k+1)
2025      {
2026        setCharacteristic (p);
2027        Variable vBuf= rootOf (gf_mipo);
2028        A= GF2FalphaRep (A, vBuf);
2029        setCharacteristic (p, extensionDeg, 'Z');
2030        ExtensionInfo info= ExtensionInfo (extension);
2031        factors= multiFactorize (A.mapinto(), info);
2032      }
2033      else // not able to pass to another GF, pass to F_p(\alpha)
2034      {
2035        setCharacteristic (p);
2036        Variable vBuf= rootOf (gf_mipo);
2037        A= GF2FalphaRep (A, vBuf);
2038        Variable v= chooseExtension (A, beta);
2039        ExtensionInfo info= ExtensionInfo (v, extension);
2040        factors= multiFactorize (A, info);
2041      }
2042    }
2043    else // need factorization over GF (p^k)
2044    {
2045      if (pow ((double) p, (double) 2*extensionDeg) < (1<<16))
2046      // pass to GF(p^2k)
2047      {
2048        setCharacteristic (p, 2*extensionDeg, 'Z');
2049        ExtensionInfo info= ExtensionInfo (k, cGFName, extension);
2050        factors= multiFactorize (GFMapUp (A, extensionDeg), info);
2051        setCharacteristic (p, extensionDeg, cGFName);
2052      }
2053      else // not able to pass to GF (p^2k), pass to F_p (\alpha)
2054      {
2055        setCharacteristic (p);
2056        Variable v1= rootOf (gf_mipo);
2057        A= GF2FalphaRep (A, v1);
2058        Variable v2= chooseExtension (A, v1);
2059        CanonicalForm primElem, imPrimElem;
2060        bool primFail= false;
2061        Variable vBuf;
2062        primElem= primitiveElement (v1, vBuf, primFail);
2063        if (primFail)
2064        {
2065          ; //ERROR
2066        }
2067        else
2068        {
2069          imPrimElem= mapPrimElem (primElem, vBuf, v2);
2070        }
2071        CFList source, dest;
2072        CanonicalForm bufA= mapUp (A, alpha, beta, primElem, imPrimElem,
2073                                     source, dest);
2074        ExtensionInfo info= ExtensionInfo (v2, v1, imPrimElem, primElem);
2075        factors= multiFactorize (bufA, info);
2076        setCharacteristic (p, k, cGFName);
2077        for (CFListIterator i= factors; i.hasItem(); i++)
2078          i.getItem()= Falpha2GFRep (i.getItem());
2079      }
2080    }
2081    return factors;
2082  }
2083}
2084
2085#endif
2086/* HAVE_NTL */
2087
Note: See TracBrowser for help on using the repository browser.