source: git/factory/facFqFactorize.cc @ b1dfaf

spielwiese
Last change on this file since b1dfaf was fd5b3a, checked in by Martin Lee <martinlee84@…>, 14 years ago
fixed logarithm git-svn-id: file:///usr/local/Singular/svn/trunk@12929 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 54.0 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  double logarithm= (double) ilog2 (totaldegree (A));
508  logarithm /= log2exp;
509  logarithm= ceil (logarithm);
510  if (factorNums < (int) logarithm)
511    factorNums= (int) logarithm;
512  for (int i= 0; i < factorNums; i++)
513  {
514    if (i == 0)
515      Aeval= A (bivarEval, y);
516    else if (i > 0 && contentAx.inCoeffDomain())
517    {
518      Aeval= A;
519      evaluation= evalPoint (A, Aeval, alpha, list, GF, fail);
520    }
521
522    if (fail && (i != 0))
523      break;
524
525    // univariate factorization
526    uniFactors= uniFactorizer (Aeval, alpha, GF);
527
528    if (uniFactors.length() == 1)
529    {
530      append (factors, contentAyFactors);
531      if (contentAyFactors.isEmpty())
532        factors.append (F/lc(F));
533      else
534      {
535        A= A (y - bivarEval, y);
536        A= A/lc (A);
537        if (!LC(A, 1).isOne())
538        {
539          CanonicalForm s,t;
540          (void) extgcd (LC (A, 1), power (y, liftBound), s, t);
541          A *= s;
542          A= mod (A, power (y, liftBound));
543        }
544        factors.append (A);
545      }
546      decompress (factors, N);
547      bivarEval= N (bivarEval);
548      return factors;
549    }
550
551    // degree analysis
552    degs= DegreePattern (uniFactors);
553
554    if (i == 0)
555    {
556      bufAeval= Aeval;
557      bufEvaluation= bivarEval;
558      bufUniFactors= uniFactors;
559      bufDegs= degs;
560      if (!contentAx.inCoeffDomain())
561        break;
562    }
563    else
564    {
565      bufDegs.intersect (degs);
566      if (uniFactors.length() < bufUniFactors.length())
567      {
568        bufUniFactors= uniFactors;
569        bufAeval= Aeval;
570        bufEvaluation= evaluation;
571      }
572    }
573
574    if (bufDegs.getLength() == 1)
575    {
576      append (factors, contentAyFactors);
577      if (contentAyFactors.isEmpty())
578        factors.append (F/lc(F));
579      else
580      {
581        A= A (y - bivarEval, y);
582        A= A/lc (A);
583        if (!LC(A, 1).isOne())
584        {
585          CanonicalForm s,t;
586          (void) extgcd (LC (A, 1), power (y, liftBound), s, t);
587          A *= s;
588          A= mod (A, power (y, liftBound));
589        }
590        factors.append (A);
591      }
592      decompress (factors, N);
593      bivarEval= N (bivarEval);
594      return factors;
595    }
596    list.append (evaluation);
597  }
598
599  bivarEval= y - bufEvaluation;
600  A= A (y + bufEvaluation, y);
601  bufUniFactors.insert (LC (A, x));
602
603  // Hensel lifting
604  CFList diophant;
605  CFMatrix M= CFMatrix (liftBound, bufUniFactors.length() - 1);
606  CFArray Pi;
607  bool earlySuccess= false;
608  int newLiftBound= 0;
609  CFList earlyFactors;
610  int smallFactorDeg= 11; //tunable parameter
611  if (smallFactorDeg >= liftBound)
612    henselLift12 (A, bufUniFactors, liftBound, Pi, diophant, M);
613  else if (smallFactorDeg >= degree (A, y) + 1)
614  {
615    henselLift12 (A, bufUniFactors, degree (A, y) + 1, Pi, diophant, M);
616    earlyFactors= earlyMonicFactorDetect (A, bufUniFactors, newLiftBound,
617                                           bufDegs, earlySuccess,
618                                           degree (A, y) + 1, liftBound);
619    if (bufDegs.getLength() > 1 && !earlySuccess)
620    {
621      if (newLiftBound > degree (A, y) + 1)
622      {
623        liftBound= newLiftBound;
624        bufUniFactors.insert (LC(A, x));
625        henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
626                            Pi, diophant, M);
627      }
628    }
629    else if (earlySuccess)
630      liftBound= newLiftBound;
631  }
632  else if (smallFactorDeg < degree (A, y) + 1)
633  {
634    henselLift12 (A, bufUniFactors, smallFactorDeg, Pi, diophant, M);
635    earlyFactors= earlyMonicFactorDetect (A, bufUniFactors, newLiftBound,
636                                           bufDegs, earlySuccess,
637                                           smallFactorDeg, liftBound);
638    if (bufDegs.getLength() > 1 && !earlySuccess)
639    {
640      bufUniFactors.insert (LC (A, x));
641      henselLiftResume12 (A, bufUniFactors, smallFactorDeg, degree (A, y) +
642                          1, Pi, diophant, M);
643      earlyFactors= earlyMonicFactorDetect (A, bufUniFactors, newLiftBound,
644                                             bufDegs, earlySuccess,
645                                             degree (A, y) + 1, liftBound);
646      if (bufDegs.getLength() > 1 && !earlySuccess)
647      {
648        if (newLiftBound > degree (A, y) + 1)
649        {
650          if (newLiftBound < newLiftBound)
651            liftBound= newLiftBound;
652          bufUniFactors.insert (LC(A, x));
653          henselLiftResume12 (A, bufUniFactors, degree (A, y) + 1, liftBound,
654                              Pi, diophant, M);
655        }
656      }
657      else if (earlySuccess)
658        liftBound= newLiftBound;
659    }
660    else if (earlySuccess)
661      liftBound= newLiftBound;
662  }
663
664  if (newLiftBound > 0)
665    liftBound= newLiftBound;
666
667  CanonicalForm MODl= power (y, liftBound);
668
669  if (bufDegs.getLength() > 1)
670    factors= monicFactorRecombi (bufUniFactors, A, MODl, bufDegs);
671
672  if (earlySuccess)
673    factors= Union (earlyFactors, factors);
674  else if (!earlySuccess && bufDegs.getLength() == 1)
675    factors= earlyFactors;
676
677  decompressAppend (factors, contentAyFactors, N);
678
679  bivarEval= N (bivarEval);
680
681  return factors;
682}
683
684CFList
685extFactorRecombination (const CFList& factors, const CanonicalForm& F,
686                        const CFList& M, const ExtensionInfo& info,
687                        const CFList& evaluation)
688{
689  Variable alpha= info.getAlpha();
690  Variable beta= info.getBeta();
691  CanonicalForm gamma= info.getGamma();
692  CanonicalForm delta= info.getDelta();
693  int k= info.getGFDegree();
694  CFList source, dest;
695  if (factors.length() == 1)
696  {
697    CanonicalForm buf= reverseShift (F, evaluation);
698    return CFList (mapDown (buf, info, source, dest));
699  }
700  if (factors.length() < 1)
701    return CFList();
702
703  int degMipoBeta;
704  if (!k && beta == Variable(1))
705    degMipoBeta= 1;
706  else if (!k && beta != Variable(1))
707    degMipoBeta= degree (getMipo (beta));
708
709  CFList T, S;
710  T= factors;
711
712  int s= 1;
713  CFList result;
714  CanonicalForm buf;
715  if (beta != Variable (1))
716    buf= mapDown (F, gamma, delta, alpha, source, dest);
717  else
718    buf= F;
719
720  CanonicalForm g, LCBuf= LC (buf, Variable (1));
721  CanonicalForm buf2;
722  int v [T.length()];
723  for (int i= 0; i < T.length(); i++)
724    v[i]= 0;
725  bool noSubset= false;
726  CFArray TT;
727  int subsetDeg;
728  TT= copy (factors);
729  bool recombination= false;
730  while (T.length() >= 2*s)
731  {
732    while (noSubset == false)
733    {
734      if (T.length() == s)
735      {
736        if (recombination)
737        {
738          T.insert (LCBuf);
739          g= prodMod (T, M);
740          T.removeFirst();
741          result.append (g/myContent (g));
742          g= reverseShift (g, evaluation);
743          g /= Lc (g);
744          appendTestMapDown (result, g, info, source, dest);
745          return result;
746        }
747        else
748        {
749          buf= reverseShift (buf, evaluation);
750          return CFList (buf);
751        }
752        return result;
753      }
754
755      S= subset (v, s, TT, noSubset);
756      if (noSubset) break;
757
758      S.insert (LCBuf);
759      g= prodMod (S, M);
760      S.removeFirst();
761      g /= myContent (g);
762      if (fdivides (g, buf))
763      {
764        buf2= reverseShift (g, evaluation);
765        buf2 /= Lc (buf2);
766        if (!k && beta == Variable (1))
767        {
768          if (degree (buf2, alpha) < degMipoBeta)
769          {
770            appendTestMapDown (result, buf2, info, source, dest);
771            buf /= g;
772            LCBuf= LC (buf, Variable (1));
773            recombination= true;
774          }
775        }
776        else
777        {
778          if (!isInExtension (buf2, delta, k))
779          {
780            appendTestMapDown (result, buf2, info, source, dest);
781            buf /= g;
782            LCBuf= LC (buf, Variable (1));
783            recombination= true;
784          }
785        }
786        T= Difference (T, S);
787
788        if (T.length() < 2*s || T.length() == s)
789        {
790          buf= reverseShift (buf, evaluation);
791          buf /= Lc (buf);
792          appendTestMapDown (result, buf, info, source, dest);
793          return result;
794        }
795        TT= copy (T);
796        indexUpdate (v, s, T.length(), noSubset);
797        if (noSubset) break;
798      }
799    }
800    s++;
801    if (T.length() < 2*s || T.length() == s)
802    {
803      buf= reverseShift (buf, evaluation);
804      appendTestMapDown (result, buf, info, source, dest);
805      return result;
806    }
807    int v [T.length()];
808    for (int i= 0; i < T.length(); i++)
809      v[i]= 0;
810    noSubset= false;
811  }
812  if (T.length() < 2*s)
813  {
814    buf= reverseShift (F, evaluation);
815    appendMapDown (result, buf, info, source, dest);
816  }
817
818  return result;
819}
820
821CFList
822factorRecombination (const CanonicalForm& F, const CFList& factors,
823                     const CFList& M)
824{
825  if (factors.length() == 1)
826    return CFList(F);
827  if (factors.length() < 1)
828    return CFList();
829
830  CFList T, S;
831
832  T= factors;
833
834  int s= 1;
835  CFList result;
836  CanonicalForm LCBuf= LC (F, Variable (1));
837  CanonicalForm g, buf= F;
838  int v [T.length()];
839  for (int i= 0; i < T.length(); i++)
840    v[i]= 0;
841  bool noSubset= false;
842  CFArray TT;
843  int subsetDeg;
844  TT= copy (factors);
845  Variable y= F.level() - 1;
846  bool recombination= false;
847  CanonicalForm h;
848  while (T.length() >= 2*s)
849  {
850    while (noSubset == false)
851    {
852      if (T.length() == s)
853      {
854        if (recombination)
855        {
856          T.insert (LC (buf));
857          g= prodMod (T, M);
858          result.append (g/myContent (g));
859          return result;
860        }
861        else
862          return CFList (F);
863      }
864      S= subset (v, s, TT, noSubset);
865      if (noSubset) break;
866      S.insert (LCBuf);
867      g= prodMod (S, M);
868      S.removeFirst();
869      g /= myContent (g);
870      if (fdivides (g, buf))
871      {
872        recombination= true;
873        result.append (g);
874        buf /= g;
875        LCBuf= LC (buf, Variable(1));
876        T= Difference (T, S);
877        if (T.length() < 2*s || T.length() == s)
878        {
879          result.append (buf);
880          return result;
881        }
882        TT= copy (T);
883        indexUpdate (v, s, T.length(), noSubset);
884        if (noSubset) break;
885      }
886    }
887    s++;
888    if (T.length() < 2*s || T.length() == s)
889    {
890      result.append (buf);
891      return result;
892    }
893    int v [T.length()];
894    for (int i= 0; i < T.length(); i++)
895      v[i]= 0;
896    noSubset= false;
897  }
898  if (T.length() < 2*s)
899    result.append (F);
900
901  return result;
902}
903
904int
905liftBoundAdaption (const CanonicalForm& F, const CFList& factors, bool&
906                   success, const int deg, const CFList& MOD, const int bound)
907{
908  int adaptedLiftBound= 0;
909  CanonicalForm buf= F;
910  Variable y= F.mvar();
911  CanonicalForm LCBuf= LC (buf, Variable (1));
912  CanonicalForm g;
913  CFList M= MOD;
914  M.append (power (y, deg));
915  int d= bound;
916  int e= 0;
917  int nBuf;
918  for (CFListIterator i= factors; i.hasItem(); i++)
919  {
920    g= mulMod (i.getItem(), LCBuf, M);
921    g /= myContent (g);
922    if (fdivides (g, buf))
923    {
924      nBuf= degree (g, y) + degree (LC (g, 1), y);
925      d -= nBuf;
926      e= tmax (e, nBuf);
927      buf /= g;
928      LCBuf= LC (buf, Variable (1));
929    }
930  }
931  adaptedLiftBound= d;
932
933  if (adaptedLiftBound < deg)
934  {
935    if (adaptedLiftBound < degree (F) + 1)
936    {
937      if (d == 1)
938      {
939        if (e + 1 > deg)
940        {
941          adaptedLiftBound= deg;
942          success= false;
943        }
944        else
945        {
946          success= true;
947          if (e + 1 < degree (F) + 1)
948            adaptedLiftBound= deg;
949          else
950            adaptedLiftBound= e + 1;
951        }
952      }
953      else
954      {
955        success= true;
956        adaptedLiftBound= deg;
957      }
958    }
959    else
960    {
961      success= true;
962    }
963  }
964  return adaptedLiftBound;
965}
966
967int
968extLiftBoundAdaption (const CanonicalForm& F, const CFList& factors, bool&
969                      success, const ExtensionInfo& info, const CFList& eval,
970                      const int deg, const CFList& MOD, const int bound)
971{
972  Variable alpha= info.getAlpha();
973  Variable beta= info.getBeta();
974  CanonicalForm gamma= info.getGamma();
975  CanonicalForm delta= info.getDelta();
976  int k= info.getGFDegree();
977  int adaptedLiftBound= 0;
978  CanonicalForm buf= F;
979  Variable y= F.mvar();
980  CanonicalForm LCBuf= LC (buf, Variable (1));
981  CanonicalForm g, gg;
982  CFList M= MOD;
983  M.append (power (y, deg));
984  adaptedLiftBound= 0;
985  int d= bound;
986  int e= 0;
987  int nBuf;
988  int degMipoBeta;
989  if (!k && beta == Variable(1))
990    degMipoBeta= 1;
991  else if (!k && beta != Variable(1))
992    degMipoBeta= degree (getMipo (beta));
993
994  for (CFListIterator i= factors; i.hasItem(); i++)
995  {
996    g= mulMod (i.getItem(), LCBuf, M);
997    g /= myContent (g);
998    if (fdivides (g, buf))
999    {
1000      gg= reverseShift (g, eval);
1001      gg /= Lc (gg);
1002      if (!k && beta == Variable (1))
1003      {
1004        if (degree (gg, alpha) < degMipoBeta)
1005        {
1006          buf /= g;
1007          nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1008          d -= nBuf;
1009          e= tmax (e, nBuf);
1010          LCBuf= LC (buf, Variable (1));
1011        }
1012      }
1013      else
1014      {
1015        if (!isInExtension (gg, delta, k))
1016        {
1017          buf /= g;
1018          nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1019          d -= nBuf;
1020          e= tmax (e, nBuf);
1021          LCBuf= LC (buf, Variable (1));
1022        }
1023      }
1024    }
1025  }
1026  adaptedLiftBound= d;
1027
1028  if (adaptedLiftBound < deg)
1029  {
1030    if (adaptedLiftBound < degree (F) + 1)
1031    {
1032      if (d == 1)
1033      {
1034        if (e + 1 > deg)
1035        {
1036          adaptedLiftBound= deg;
1037          success= false;
1038        }
1039        else
1040        {
1041          success= true;
1042          if (e + 1 < degree (F) + 1)
1043            adaptedLiftBound= deg;
1044          else
1045            adaptedLiftBound= e + 1;
1046        }
1047      }
1048      else
1049      {
1050        success= true;
1051        adaptedLiftBound= deg;
1052      }
1053    }
1054    else
1055    {
1056      success= true;
1057    }
1058  }
1059
1060  return adaptedLiftBound;
1061}
1062
1063CFList
1064earlyFactorDetect (CanonicalForm& F, CFList& factors, int& adaptedLiftBound,
1065                   bool& success, const int deg, const CFList& MOD,
1066                   const int bound)
1067{
1068  CFList result;
1069  CFList T= factors;
1070  CanonicalForm buf= F;
1071  Variable y= F.mvar();
1072  CanonicalForm LCBuf= LC (buf, Variable (1));
1073  CanonicalForm g;
1074  CFList M= MOD;
1075  M.append (power (y, deg));
1076  adaptedLiftBound= 0;
1077  int d= bound;
1078  int e= 0;
1079  int nBuf;
1080  for (CFListIterator i= factors; i.hasItem(); i++)
1081  {
1082    g= mulMod (i.getItem(), LCBuf, M);
1083    g /= myContent (g);
1084    if (fdivides (g, buf))
1085    {
1086      result.append (g);
1087      nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1088      d -= nBuf;
1089      e= tmax (e, nBuf);
1090      buf /= g;
1091      LCBuf= LC (buf, Variable (1));
1092      T= Difference (T, CFList (i.getItem()));
1093    }
1094  }
1095  adaptedLiftBound= d;
1096
1097  if (adaptedLiftBound < deg)
1098  {
1099    if (adaptedLiftBound < degree (F) + 1)
1100    {
1101      if (d == 1)
1102        adaptedLiftBound= tmin (e + 1, deg);
1103      else
1104        adaptedLiftBound= deg;
1105    }
1106    factors= T;
1107    F= buf;
1108    success= true;
1109  }
1110  return result;
1111}
1112
1113CFList
1114extEarlyFactorDetect (CanonicalForm& F, CFList& factors, int& adaptedLiftBound,
1115                      bool& success, const ExtensionInfo& info, const CFList&
1116                      eval, const int deg, const CFList& MOD, const int bound)
1117{
1118  Variable alpha= info.getAlpha();
1119  Variable beta= info.getBeta();
1120  CanonicalForm gamma= info.getGamma();
1121  CanonicalForm delta= info.getDelta();
1122  int k= info.getGFDegree();
1123  CFList result;
1124  CFList T= factors;
1125  CanonicalForm buf= F;
1126  Variable y= F.mvar();
1127  CanonicalForm LCBuf= LC (buf, Variable (1));
1128  CanonicalForm g, gg;
1129  CFList M= MOD;
1130  M.append (power (y, deg));
1131  adaptedLiftBound= 0;
1132  int d= bound;
1133  int e= 0;
1134  int nBuf;
1135  CFList source, dest;
1136
1137  int degMipoBeta;
1138  if (!k && beta == Variable(1))
1139    degMipoBeta= 1;
1140  else if (!k && beta != Variable(1))
1141    degMipoBeta= degree (getMipo (beta));
1142
1143  for (CFListIterator i= factors; i.hasItem(); i++)
1144  {
1145    g= mulMod (i.getItem(), LCBuf, M);
1146    g /= myContent (g);
1147    if (fdivides (g, buf))
1148    {
1149      gg= reverseShift (g, eval);
1150      gg /= Lc (gg);
1151          if (!k && beta == Variable (1))
1152          {
1153            if (degree (gg, alpha) < degMipoBeta)
1154            {
1155              appendTestMapDown (result, gg, info, source, dest);
1156              buf /= g;
1157              nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1158              d -= nBuf;
1159              e= tmax (e, nBuf);
1160              LCBuf= LC (buf, Variable (1));
1161            }
1162          }
1163          else
1164          {
1165            if (!isInExtension (gg, delta, k))
1166            {
1167              appendTestMapDown (result, gg, info, source, dest);
1168              buf /= g;
1169              nBuf= degree (g, y) + degree (LC (g, Variable (1)), y);
1170              d -= nBuf;
1171              e= tmax (e, nBuf);
1172              LCBuf= LC (buf, Variable (1));
1173            }
1174          }
1175      T= Difference (T, CFList (i.getItem()));
1176    }
1177  }
1178  adaptedLiftBound= d;
1179
1180  if (adaptedLiftBound < deg)
1181  {
1182    if (adaptedLiftBound < degree (F) + 1)
1183    {
1184      if (d == 1)
1185        adaptedLiftBound= tmin (e + 1, deg);
1186      else
1187        adaptedLiftBound= deg;
1188    }
1189    success= true;
1190    factors= T;
1191    F= buf;
1192  }
1193  return result;
1194}
1195
1196CFList
1197evalPoints (const CanonicalForm& F, CFList & eval, const Variable& alpha,
1198            CFList& list, const bool& GF, bool& fail)
1199{
1200  int k= F.level() - 1;
1201  Variable x= Variable (1);
1202  CFList result;
1203  FFRandom genFF;
1204  GFRandom genGF;
1205  int p= getCharacteristic ();
1206  double bound;
1207  if (alpha != Variable (1))
1208  {
1209    bound= pow ((double) p, (double) degree (getMipo(alpha)));
1210    bound= pow ((double) bound, (double) k);
1211  }
1212  else if (GF)
1213  {
1214    bound= pow ((double) p, (double) getGFDegree());
1215    bound= pow ((double) bound, (double) k);
1216  }
1217  else
1218    bound= pow ((double) p, (double) k);
1219
1220  CanonicalForm random;
1221  CanonicalForm deriv_x, gcd_deriv;
1222  do
1223  {
1224    random= 0;
1225    // possible overflow if list.length() does not fit into a int
1226    if (list.length() >= bound)
1227    {
1228      fail= true;
1229      break;
1230    }
1231    for (int i= 0; i < k; i++)
1232    {
1233      if (GF)
1234      {
1235        result.append (genGF.generate());
1236        random += result.getLast()*power (x, i);
1237      }
1238      else if (alpha != Variable(1))
1239      {
1240        AlgExtRandomF genAlgExt (alpha);
1241        result.append (genAlgExt.generate());
1242        random += result.getLast()*power (x, i);
1243      }
1244      else
1245      {
1246        result.append (genFF.generate());
1247        random += result.getLast()*power (x, i);
1248      }
1249    }
1250    if (find (list, random))
1251    {
1252      result= CFList();
1253      continue;
1254    }
1255    int l= F.level();
1256    eval.insert (F);
1257    for (CFListIterator i= result; i.hasItem(); i++, l--)
1258      eval.insert (eval.getFirst()(i.getItem(), l));
1259
1260    if (degree (eval.getFirst()) != degree (F, 1))
1261    {
1262      if (!find (list, random))
1263        list.append (random);
1264      result= CFList();
1265      eval= CFList();
1266      continue;
1267    }
1268
1269    deriv_x= deriv (eval.getFirst(), x);
1270    gcd_deriv= gcd (eval.getFirst(), deriv_x);
1271    if (degree (gcd_deriv) > 0)
1272    {
1273      if (!find (list, random))
1274        list.append (random);
1275      result= CFList();
1276      eval= CFList();
1277      continue;
1278    }
1279    CFListIterator i= eval;
1280    i++;
1281    CanonicalForm contentx= content (i.getItem(), x);
1282    if (degree (contentx) > 0)
1283    {
1284      if (!find (list, random))
1285        list.append (random);
1286      result= CFList();
1287      eval= CFList();
1288      continue;
1289    }
1290
1291    if (list.length() >= bound)
1292    {
1293      fail= true;
1294      break;
1295    }
1296  } while (find (list, random));
1297
1298  if (!eval.isEmpty())
1299    eval.removeFirst();
1300
1301  return result;
1302}
1303
1304static inline
1305int newMainVariableSearch (CanonicalForm& A, CFList& Aeval, CFList&
1306                           evaluation, const Variable& alpha, const int lev)
1307{
1308  Variable x= Variable (1);
1309  CanonicalForm derivI, buf;
1310  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1311  int swapLevel= 0;
1312  CFList list;
1313  bool fail= false;
1314  buf= A;
1315  Aeval= CFList();
1316  evaluation= CFList();
1317  for (int i= lev; i <= A.level(); i++)
1318  {
1319    derivI= deriv (buf, Variable (i));
1320    if (!derivI.isZero())
1321    {
1322      buf= swapvar (buf, x, Variable (i));
1323      Aeval= CFList();
1324      evaluation= CFList();
1325      fail= false;
1326      evaluation= evalPoints (buf, Aeval, alpha, list, GF, fail);
1327      if (!fail)
1328      {
1329        A= buf;
1330        swapLevel= i;
1331        break;
1332      }
1333      else
1334        buf= A;
1335    }
1336  }
1337  return swapLevel;
1338}
1339
1340static inline
1341CanonicalForm lcmContent (const CanonicalForm& A, CFList& contentAi)
1342{
1343  int i= A.level();
1344  contentAi.append (myContent (A, i));
1345  contentAi.append (myContent (A, i - 1));
1346  CanonicalForm result= myLcm (contentAi.getFirst(), contentAi.getLast());
1347  for (i= i - 2; i > 0; i--)
1348  {
1349    contentAi.append (content (A, i));
1350    result= myLcm (result, contentAi.getLast());
1351  }
1352  return result;
1353}
1354
1355CFList
1356henselLiftAndEarly (CanonicalForm& A, CFList& MOD, int*& liftBounds, bool&
1357                    earlySuccess, CFList& earlyFactors, const CFList& Aeval,
1358                    const CFList& biFactors, const CFList& evaluation,
1359                    const ExtensionInfo& info)
1360{
1361  bool extension= info.isInExtension();
1362  CFList bufFactors= biFactors;
1363  bufFactors.insert (LC (Aeval.getFirst(), 1));
1364
1365  sortList (bufFactors, Variable (1));
1366
1367  CFList diophant;
1368  CFArray Pi;
1369  int smallFactorDeg= 11; //tunable parameter
1370  CFList result;
1371  int newLiftBound= 0;
1372  int adaptedLiftBound= 0;
1373  int liftBound= liftBounds[1];
1374
1375  earlySuccess= false;
1376  bool earlyReconst= false;
1377  CFList earlyReconstFactors;
1378  CFListIterator j= Aeval;
1379  j++;
1380  CanonicalForm buf= j.getItem();
1381  CFMatrix Mat= CFMatrix (liftBound, bufFactors.length() - 1);
1382  MOD= CFList (power (Variable (2), liftBounds[0]));
1383  if (smallFactorDeg >= liftBound)
1384  {
1385    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
1386  }
1387  else if (smallFactorDeg >= degree (buf) + 1)
1388  {
1389    liftBounds[1]= degree (buf) + 1;
1390    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
1391    if (Aeval.length() == 2)
1392    {
1393      if (!extension)
1394        earlyFactors= earlyFactorDetect
1395                       (buf, result, adaptedLiftBound, earlySuccess,
1396                        degree (buf) + 1, MOD, liftBound);
1397      else
1398        earlyFactors= extEarlyFactorDetect
1399                       (buf, result, adaptedLiftBound, earlySuccess,
1400                        info, evaluation, degree
1401                        (buf) + 1, MOD, liftBound);
1402    }
1403    else
1404    {
1405      if (!extension)
1406        adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
1407                                             degree (buf) + 1, MOD, liftBound);
1408      else
1409        adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess, info,
1410                                                evaluation, degree (buf) + 1,
1411                                                MOD, liftBound);
1412    }
1413    if (!earlySuccess)
1414    {
1415      result.insert (LC (buf, 1));
1416      liftBounds[1]= adaptedLiftBound;
1417      liftBound= adaptedLiftBound;
1418      henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1419                        Pi, diophant, Mat, MOD);
1420    }
1421    else
1422      liftBounds[1]= adaptedLiftBound;
1423  }
1424  else if (smallFactorDeg < degree (buf) + 1)
1425  {
1426    liftBounds[1]= smallFactorDeg;
1427    result= henselLift23 (Aeval, bufFactors, liftBounds, diophant, Pi, Mat);
1428    if (Aeval.length() == 2)
1429    {
1430      if (!extension)
1431        earlyFactors= earlyFactorDetect (buf, result, adaptedLiftBound,
1432                                         earlySuccess, smallFactorDeg, MOD,
1433                                         liftBound);
1434      else
1435        earlyFactors= extEarlyFactorDetect (buf, result, adaptedLiftBound,
1436                                            earlySuccess, info, evaluation,
1437                                            smallFactorDeg, MOD, liftBound);
1438    }
1439    else
1440    {
1441      if (!extension)
1442        adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
1443                                             smallFactorDeg, MOD, liftBound);
1444      else
1445        adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess, info,
1446                                                evaluation, smallFactorDeg, MOD,
1447                                                liftBound);
1448    }
1449
1450    if (!earlySuccess)
1451    {
1452      result.insert (LC (buf, 1));
1453      henselLiftResume (buf, result, smallFactorDeg, degree (buf) + 1,
1454                        Pi, diophant, Mat, MOD);
1455      if (Aeval.length() == 2)
1456      {
1457         if (!extension)
1458           earlyFactors= earlyFactorDetect (buf, result, adaptedLiftBound,
1459                                            earlySuccess, degree (buf) + 1,
1460                                            MOD, liftBound);
1461         else
1462           earlyFactors= extEarlyFactorDetect (buf, result, adaptedLiftBound,
1463                                               earlySuccess, info, evaluation,
1464                                               degree (buf) + 1, MOD,
1465                                               liftBound);
1466      }
1467      else
1468      {
1469        if (!extension)
1470          adaptedLiftBound= liftBoundAdaption (buf, result, earlySuccess,
1471                                               degree (buf) + 1, MOD,liftBound);
1472        else
1473          adaptedLiftBound= extLiftBoundAdaption (buf, result, earlySuccess,
1474                                                  info, evaluation,
1475                                                  degree (buf) + 1, MOD,
1476                                                  liftBound);
1477      }
1478      if (!earlySuccess)
1479      {
1480        result.insert (LC (buf, 1));
1481        liftBounds[1]= adaptedLiftBound;
1482        liftBound= adaptedLiftBound;
1483        henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1484                          Pi, diophant, Mat, MOD);
1485      }
1486      else
1487        liftBounds[1]= adaptedLiftBound;
1488    }
1489    else
1490      liftBounds[1]= adaptedLiftBound;
1491  }
1492
1493  MOD.append (power (Variable (3), liftBounds[1]));
1494
1495  if (Aeval.length() > 2)
1496  {
1497    CFListIterator j= Aeval;
1498    j++;
1499    CFList bufEval;
1500    bufEval.append (j.getItem());
1501    j++;
1502    int liftBoundsLength= Aeval.getLast().level() - 1;
1503    for (int i= 2; i <= liftBoundsLength && j.hasItem(); i++, j++)
1504    {
1505      earlySuccess= false;
1506      result.insert (LC (bufEval.getFirst(), 1));
1507      bufEval.append (j.getItem());
1508      liftBound= liftBounds[i];
1509      Mat= CFMatrix (liftBounds[i], result.length() - 1);
1510
1511      buf= j.getItem();
1512      if (smallFactorDeg >= liftBound)
1513        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
1514                            liftBounds[i -  1], liftBounds[i]);
1515      else if (smallFactorDeg >= degree (buf) + 1)
1516      {
1517        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
1518                            liftBounds[i -  1], degree (buf) + 1);
1519
1520        if (Aeval.length() == i + 1)
1521        {
1522          if (!extension)
1523            earlyFactors= earlyFactorDetect
1524                           (buf, result, adaptedLiftBound, earlySuccess,
1525                            degree (buf) + 1, MOD, liftBound);
1526          else
1527            earlyFactors= extEarlyFactorDetect
1528                           (buf, result, adaptedLiftBound, earlySuccess,
1529                            info, evaluation, degree (buf) + 1, MOD, liftBound);
1530        }
1531        else
1532        {
1533          if (!extension)
1534            adaptedLiftBound= liftBoundAdaption
1535                                (buf, result, earlySuccess, degree (buf)
1536                                 + 1,  MOD, liftBound);
1537          else
1538            adaptedLiftBound= extLiftBoundAdaption
1539                                (buf, result, earlySuccess, info, evaluation,
1540                                 degree (buf) + 1, MOD, liftBound);
1541        }
1542
1543        if (!earlySuccess)
1544        {
1545          result.insert (LC (buf, 1));
1546          liftBounds[i]= adaptedLiftBound;
1547          liftBound= adaptedLiftBound;
1548          henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1549                            Pi, diophant, Mat, MOD);
1550        }
1551        else
1552        {
1553          liftBounds[i]= adaptedLiftBound;
1554        }
1555      }
1556      else if (smallFactorDeg < degree (buf) + 1)
1557      {
1558        result= henselLift (bufEval, result, MOD, diophant, Pi, Mat,
1559                            liftBounds[i -  1], smallFactorDeg);
1560
1561        if (Aeval.length() == i + 1)
1562        {
1563          if (!extension)
1564            earlyFactors= earlyFactorDetect
1565                           (buf, result, adaptedLiftBound, earlySuccess,
1566                            smallFactorDeg, MOD, liftBound);
1567          else
1568            earlyFactors= extEarlyFactorDetect
1569                           (buf, result, adaptedLiftBound, earlySuccess,
1570                            info, evaluation, smallFactorDeg, MOD, liftBound);
1571        }
1572        else
1573        {
1574          if (!extension)
1575            adaptedLiftBound= liftBoundAdaption
1576                                (buf, result, earlySuccess,
1577                                 smallFactorDeg, MOD, liftBound);
1578          else
1579            adaptedLiftBound= extLiftBoundAdaption
1580                                (buf, result, earlySuccess, info, evaluation,
1581                                 smallFactorDeg, MOD, liftBound);
1582        }
1583
1584        if (!earlySuccess)
1585        {
1586          result.insert (LC (buf, 1));
1587          henselLiftResume (buf, result, smallFactorDeg,
1588                            degree (buf) + 1, Pi, diophant, Mat, MOD);
1589          if (Aeval.length() == i + 1)
1590          {
1591            if (!extension)
1592              earlyFactors= earlyFactorDetect
1593                             (buf, result, adaptedLiftBound, earlySuccess,
1594                              degree (buf) +  1,  MOD, liftBound);
1595            else
1596              earlyFactors= extEarlyFactorDetect
1597                             (buf, result, adaptedLiftBound, earlySuccess,
1598                              info, evaluation, degree (buf) + 1, MOD,
1599                              liftBound);
1600          }
1601          else
1602          {
1603            if (!extension)
1604              adaptedLiftBound= liftBoundAdaption
1605                                  (buf, result, earlySuccess, degree
1606                                   (buf) +  1,  MOD, liftBound);
1607            else
1608              adaptedLiftBound= extLiftBoundAdaption
1609                                  (buf, result, earlySuccess, info, evaluation,
1610                                   degree (buf) + 1,  MOD, liftBound);
1611          }
1612
1613          if (!earlySuccess)
1614          {
1615            result.insert (LC (buf, 1));
1616            liftBounds[i]= adaptedLiftBound;
1617            liftBound= adaptedLiftBound;
1618            henselLiftResume (buf, result, degree (buf) + 1, liftBound,
1619                              Pi, diophant, Mat, MOD);
1620          }
1621          else
1622            liftBounds[i]= adaptedLiftBound;
1623        }
1624        else
1625          liftBounds[i]= adaptedLiftBound;
1626      }
1627      MOD.append (power (Variable (i + 2), liftBounds[i]));
1628      bufEval.removeFirst();
1629    }
1630    bufFactors= result;
1631  }
1632  else
1633    bufFactors= result;
1634
1635  if (earlySuccess)
1636    A= buf;
1637  return result;
1638}
1639
1640CFList
1641extFactorize (const CanonicalForm& F, const ExtensionInfo& info);
1642
1643CFList
1644multiFactorize (const CanonicalForm& F, const ExtensionInfo& info)
1645{
1646
1647  if (F.inCoeffDomain())
1648    return CFList (F);
1649
1650  // compress and find main Variable
1651  CFMap N;
1652  CanonicalForm A= myCompress (F, N);
1653
1654  A /= Lc (A); // make monic
1655
1656  Variable alpha= info.getAlpha();
1657  Variable beta= info.getBeta();
1658  CanonicalForm gamma= info.getGamma();
1659  CanonicalForm delta= info.getDelta();
1660  int k= info.getGFDegree();
1661  bool extension= info.isInExtension();
1662  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1663  //univariate case
1664  if (F.isUnivariate())
1665  {
1666    if (extension == false)
1667      return uniFactorizer (F, alpha, GF);
1668    else
1669    {
1670      CFList source, dest;
1671      A= mapDown (F, info, source, dest);
1672      return uniFactorizer (A, beta, GF);
1673    }
1674  }
1675
1676  //bivariate case
1677  if (A.level() == 2)
1678  {
1679    CFList buf= biFactorize (F, info);
1680    return buf;
1681  }
1682
1683  Variable x= Variable (1);
1684  Variable y= Variable (2);
1685
1686  // remove content
1687  CFList contentAi;
1688  CanonicalForm lcmCont= lcmContent (A, contentAi);
1689  A /= lcmCont;
1690
1691  // trivial after content removal
1692  CFList contentAFactors;
1693  if (A.inCoeffDomain())
1694  {
1695    for (CFListIterator i= contentAi; i.hasItem(); i++)
1696    {
1697      if (i.getItem().inCoeffDomain())
1698        continue;
1699      else
1700      {
1701        lcmCont /= i.getItem();
1702        contentAFactors=
1703        Union (multiFactorize (lcmCont, info),
1704               multiFactorize (i.getItem(), info));
1705        break;
1706      }
1707    }
1708    decompress (contentAFactors, N);
1709    normalize (contentAFactors);
1710    return contentAFactors;
1711  }
1712
1713  // factorize content
1714  contentAFactors= multiFactorize (lcmCont, info);
1715
1716  // univariate after content removal
1717  CFList factors;
1718  if (A.isUnivariate ())
1719  {
1720    factors= uniFactorizer (A, alpha, GF);
1721    append (factors, contentAFactors);
1722    decompress (factors, N);
1723    return factors;
1724  }
1725
1726  // check main variable
1727  int swapLevel= 0;
1728  CanonicalForm derivZ;
1729  CanonicalForm gcdDerivZ;
1730  CanonicalForm bufA= A;
1731  Variable z;
1732  for (int i= 1; i <= A.level(); i++)
1733  {
1734    z= Variable (i);
1735    derivZ= deriv (bufA, z);
1736    if (derivZ.isZero())
1737    {
1738      if (i == 1)
1739        swapLevel= 1;
1740      else
1741        continue;
1742    }
1743    else
1744    {
1745      if (swapLevel == 1)
1746      {
1747        swapLevel= i;
1748        A= swapvar (A, x, z);
1749      }
1750      if (GF)
1751        gcdDerivZ= GCD_GF (bufA, derivZ);
1752      else if (alpha == Variable (1))
1753        gcdDerivZ= GCD_small_p (bufA, derivZ);
1754      else
1755        gcdDerivZ= GCD_Fp_extension (bufA, derivZ, alpha);
1756      if (degree (gcdDerivZ) > 0 && !derivZ.isZero())
1757      {
1758        CanonicalForm g= bufA/gcdDerivZ;
1759        CFList factorsG=
1760        Union (multiFactorize (g, info),
1761               multiFactorize (gcdDerivZ, info));
1762        appendSwapDecompress (factorsG, contentAFactors, N, swapLevel, x);
1763        normalize (factorsG);
1764        return factorsG;
1765      }
1766    }
1767  }
1768
1769
1770  CFList Aeval, list, evaluation, bufEvaluation, bufAeval;
1771  bool fail= false;
1772  int swapLevel2= 0;
1773  int level;
1774  int factorNums= 3;
1775  CanonicalForm bivarEval;
1776  CFList biFactors, bufBiFactors;
1777  CanonicalForm evalPoly;
1778  int lift, bufLift;
1779  double logarithm= (double) ilog2 (totaldegree (A));
1780  logarithm /= log2exp;
1781  logarithm= ceil (logarithm);
1782  if (factorNums < (int) logarithm)
1783    factorNums= (int) logarithm;
1784  // several bivariate factorizations
1785  for (int i= 0; i < factorNums; i++)
1786  {
1787    bufA= A;
1788    bufAeval= CFList();
1789    bufEvaluation= evalPoints (bufA, bufAeval, alpha, list, GF, fail);
1790    evalPoly= 0;
1791
1792    if (fail && (i == 0))
1793    {
1794      if (!swapLevel)
1795        level= 2;
1796      else
1797        level= swapLevel + 1;
1798
1799      swapLevel2= newMainVariableSearch (A, Aeval, evaluation, alpha, level);
1800
1801      if (!swapLevel2) // need to pass to an extension
1802      {
1803        factors= extFactorize (A, info);
1804        appendSwapDecompress (factors, contentAFactors, N, swapLevel, x);
1805        normalize (factors);
1806        return factors;
1807      }
1808      else
1809      {
1810        fail= false;
1811        bufAeval= Aeval;
1812        bufA= A;
1813        bufEvaluation= evaluation;
1814      }
1815    }
1816    else if (fail && (i > 0))
1817      break;
1818
1819    bivarEval= bufEvaluation.getLast();
1820    bufEvaluation.removeLast();
1821
1822    bufLift= degree (A, y) + 1 + degree (LC(A, x), y);
1823
1824    TIMING_START (fac_bi_factorizer);
1825    bufBiFactors= biFactorizer (bufAeval.getFirst(), alpha, bivarEval, bufLift);
1826    TIMING_END_AND_PRINT (fac_bi_factorizer,
1827                          "time for bivariate factorization: ");
1828
1829    if (bufBiFactors.length() == 1)
1830    {
1831      if (extension)
1832      {
1833        CFList source, dest;
1834        A= mapDown (A, info, source, dest);
1835      }
1836      factors.append (A);
1837      appendSwapDecompress (factors, contentAFactors, N, swapLevel,
1838                            swapLevel2, x);
1839      normalize (factors);
1840      return factors;
1841    }
1842
1843    bufEvaluation.append (-bivarEval[0]);
1844    if (i == 0)
1845    {
1846      Aeval= bufAeval;
1847      evaluation= bufEvaluation;
1848      biFactors= bufBiFactors;
1849      lift= bufLift;
1850    }
1851    else
1852    {
1853      if (bufBiFactors.length() < biFactors.length() ||
1854          ((bufLift < lift) && (bufBiFactors.length() == biFactors.length())))
1855      {
1856        Aeval= bufAeval;
1857        evaluation= bufEvaluation;
1858        biFactors= bufBiFactors;
1859        lift= bufLift;
1860      }
1861    }
1862    int k= 0;
1863    for (CFListIterator j= bufEvaluation; j.hasItem(); j++, k++)
1864      evalPoly += j.getItem()*power (x, k);
1865    list.append (evalPoly);
1866  }
1867
1868  //shifting to zero
1869  A= shift2Zero (A, Aeval, evaluation);
1870
1871  int* liftBounds;
1872  int liftBoundsLength= F.level() - 1;
1873  liftBounds= liftingBounds (A, lift);
1874
1875  CFList MOD;
1876  bool earlySuccess;
1877  CFList earlyFactors;
1878  TIMING_START (fac_hensel_lift);
1879  CFList liftedFactors= henselLiftAndEarly
1880                        (A, MOD, liftBounds, earlySuccess, earlyFactors,
1881                         Aeval, biFactors, evaluation, info);
1882  TIMING_END_AND_PRINT (fac_hensel_lift, "time for hensel lifting: ");
1883
1884  if (!extension)
1885  {
1886    TIMING_START (fac_factor_recombination);
1887    factors= factorRecombination (A, liftedFactors, MOD);
1888    TIMING_END_AND_PRINT (fac_factor_recombination,
1889                          "time for factor recombination: ");
1890  }
1891  else
1892  {
1893    TIMING_START (fac_factor_recombination);
1894    factors= extFactorRecombination (liftedFactors, A, MOD, info, evaluation);
1895    TIMING_END_AND_PRINT (fac_factor_recombination,
1896                          "time for factor recombination: ");
1897  }
1898
1899  if (earlySuccess)
1900    factors= Union (factors, earlyFactors);
1901
1902  if (!extension)
1903  {
1904    for (CFListIterator i= factors; i.hasItem(); i++)
1905    {
1906      int kk= Aeval.getLast().level();
1907      for (CFListIterator j= evaluation; j.hasItem(); j++, kk--)
1908      {
1909        if (i.getItem().level() < kk)
1910          continue;
1911        i.getItem()= i.getItem() (Variable (kk) - j.getItem(), kk);
1912      }
1913    }
1914  }
1915
1916  swap (factors, swapLevel, swapLevel2, x);
1917  append (factors, contentAFactors);
1918  decompress (factors, N);
1919  normalize (factors);
1920
1921  delete[] liftBounds;
1922
1923  return factors;
1924}
1925
1926/// multivariate factorization over an extension of the initial field
1927CFList
1928extFactorize (const CanonicalForm& F, const ExtensionInfo& info)
1929{
1930  CanonicalForm A= F;
1931
1932  Variable alpha= info.getAlpha();
1933  Variable beta= info.getBeta();
1934  int k= info.getGFDegree();
1935  char cGFName= info.getGFName();
1936  CanonicalForm delta= info.getDelta();
1937  bool GF= (CFFactory::gettype() == GaloisFieldDomain);
1938  Variable w= Variable (1);
1939
1940  CFList factors;
1941  if (!GF && alpha == w)  // we are in F_p
1942  {
1943    CFList factors;
1944    bool extension= true;
1945    int p= getCharacteristic();
1946    if (p*p < (1<<16)) // pass to GF if possible
1947    {
1948      setCharacteristic (getCharacteristic(), 2, 'Z');
1949      ExtensionInfo info= ExtensionInfo (extension);
1950      A= A.mapinto();
1951      factors= multiFactorize (A, info);
1952
1953      Variable vBuf= rootOf (gf_mipo);
1954      setCharacteristic (getCharacteristic());
1955      for (CFListIterator j= factors; j.hasItem(); j++)
1956        j.getItem()= GF2FalphaRep (j.getItem(), vBuf);
1957    }
1958    else  // not able to pass to GF, pass to F_p(\alpha)
1959    {
1960      Variable v= chooseExtension (A, beta);
1961      ExtensionInfo info= ExtensionInfo (v, extension);
1962      factors= multiFactorize (A, info);
1963    }
1964    return factors;
1965  }
1966  else if (!GF && (alpha != w)) // we are in F_p(\alpha)
1967  {
1968    if (k == 1) // need factorization over F_p
1969    {
1970      int extDeg= degree (getMipo (alpha));
1971      extDeg++;
1972      CanonicalForm mipo= randomIrredpoly (extDeg + 1, Variable (1));
1973      Variable v= rootOf (mipo);
1974      ExtensionInfo info= ExtensionInfo (v);
1975      factors= biFactorize (A, info);
1976    }
1977    else
1978    {
1979      if (beta == Variable (1))
1980      {
1981        Variable v= chooseExtension (A, alpha);
1982        CanonicalForm primElem, imPrimElem;
1983        bool primFail= false;
1984        Variable vBuf;
1985        primElem= primitiveElement (alpha, vBuf, primFail);
1986        ASSERT (!primFail, "failure in integer factorizer");
1987        if (primFail)
1988          ; //ERROR
1989        else
1990          imPrimElem= mapPrimElem (primElem, vBuf, v);
1991
1992        CFList source, dest;
1993        CanonicalForm bufA= mapUp (A, alpha, v, primElem, imPrimElem,
1994                                   source, dest);
1995        ExtensionInfo info= ExtensionInfo (v, alpha, imPrimElem, primElem);
1996        factors= biFactorize (bufA, info);
1997      }
1998      else
1999      {
2000        Variable v= chooseExtension (A, alpha);
2001        CanonicalForm primElem, imPrimElem;
2002        bool primFail= false;
2003        Variable vBuf;
2004        ASSERT (!primFail, "failure in integer factorizer");
2005        if (primFail)
2006          ; //ERROR
2007        else
2008          imPrimElem= mapPrimElem (delta, beta, v); //oder mapPrimElem (primElem, vBuf, v);
2009
2010        CFList source, dest;
2011        CanonicalForm bufA= mapDown (A, info, source, dest);
2012        source= CFList();
2013        dest= CFList();
2014        bufA= mapUp (bufA, beta, v, delta, imPrimElem, source, dest);
2015        ExtensionInfo info= ExtensionInfo (v, beta, imPrimElem, delta);
2016        factors= biFactorize (bufA, info);
2017      }
2018    }
2019    return factors;
2020  }
2021  else // we are in GF (p^k)
2022  {
2023    int p= getCharacteristic();
2024    int extensionDeg= getGFDegree();
2025    bool extension= true;
2026    if (k == 1) // need factorization over F_p
2027    {
2028      extensionDeg++;
2029      if (pow ((double) p, (double) extensionDeg) < (1<<16))
2030      // pass to GF(p^k+1)
2031      {
2032        setCharacteristic (p);
2033        Variable vBuf= rootOf (gf_mipo);
2034        A= GF2FalphaRep (A, vBuf);
2035        setCharacteristic (p, extensionDeg, 'Z');
2036        ExtensionInfo info= ExtensionInfo (extension);
2037        factors= multiFactorize (A.mapinto(), info);
2038      }
2039      else // not able to pass to another GF, pass to F_p(\alpha)
2040      {
2041        setCharacteristic (p);
2042        Variable vBuf= rootOf (gf_mipo);
2043        A= GF2FalphaRep (A, vBuf);
2044        Variable v= chooseExtension (A, beta);
2045        ExtensionInfo info= ExtensionInfo (v, extension);
2046        factors= multiFactorize (A, info);
2047      }
2048    }
2049    else // need factorization over GF (p^k)
2050    {
2051      if (pow ((double) p, (double) 2*extensionDeg) < (1<<16))
2052      // pass to GF(p^2k)
2053      {
2054        setCharacteristic (p, 2*extensionDeg, 'Z');
2055        ExtensionInfo info= ExtensionInfo (k, cGFName, extension);
2056        factors= multiFactorize (GFMapUp (A, extensionDeg), info);
2057        setCharacteristic (p, extensionDeg, cGFName);
2058      }
2059      else // not able to pass to GF (p^2k), pass to F_p (\alpha)
2060      {
2061        setCharacteristic (p);
2062        Variable v1= rootOf (gf_mipo);
2063        A= GF2FalphaRep (A, v1);
2064        Variable v2= chooseExtension (A, v1);
2065        CanonicalForm primElem, imPrimElem;
2066        bool primFail= false;
2067        Variable vBuf;
2068        primElem= primitiveElement (v1, vBuf, primFail);
2069        if (primFail)
2070        {
2071          ; //ERROR
2072        }
2073        else
2074        {
2075          imPrimElem= mapPrimElem (primElem, vBuf, v2);
2076        }
2077        CFList source, dest;
2078        CanonicalForm bufA= mapUp (A, alpha, beta, primElem, imPrimElem,
2079                                     source, dest);
2080        ExtensionInfo info= ExtensionInfo (v2, v1, imPrimElem, primElem);
2081        factors= multiFactorize (bufA, info);
2082        setCharacteristic (p, k, cGFName);
2083        for (CFListIterator i= factors; i.hasItem(); i++)
2084          i.getItem()= Falpha2GFRep (i.getItem());
2085      }
2086    }
2087    return factors;
2088  }
2089}
2090
2091#endif
2092/* HAVE_NTL */
2093
Note: See TracBrowser for help on using the repository browser.